# Approximating Robot Configuration Spaces with few Convex Sets
using Clique Covers of Visibility Graphs

Peter Werner, Alexandre Amice, Tobia Marcucci, Daniela Rus, and Russ Tedrake
All authors are with CSAIL, MIT. The corresponding author is P.Werner,
<EMAIL_ADDRESS>This work was supported in part by Amazon.com PO#
#2D-06310236, Air Force Research Laboratory FA8750-19-2-1000, Office of Naval
Research (ONR), Awards No. N00014-22-1-2121 and N00014-18-1-2830, and the
Toyota Research Institute (TRI).

###### Abstract

Many computations in robotics can be dramatically accelerated if the robot
configuration space is described as a collection of simple sets. For example,
recently developed motion planners rely on a convex decomposition of the free
space to design collision-free trajectories using fast convex optimization. In
this work, we present an efficient method for approximately covering complex
configuration spaces with a small number of polytopes. The approach constructs
a visibility graph using sampling and generates a clique cover of this graph
to find clusters of samples that have mutual line of sight. These clusters are
then inflated into large, full-dimensional, polytopes. We evaluate our method
on a variety of robotic systems and show that it consistently covers larger
portions of free configuration space, with fewer polytopes, and in a fraction
of the time compared to previous methods.

##  I Introduction

Approximating complex sets as a union of simpler sets is a common pre-
processing method for accelerating downstream computations. Familiar examples
include clustering methods for high-dimensional data in machine learning [1,
§2.5.2], approximating complex shapes via triangular meshes to facilitate
graphics rendering [2], and describing geometries as unions of convex sets for
efficient collision checking [3]. Similarly, recently developed methods for
robot motion planning rely on (conservative) decompositions of the environment
into convex sets to design smooth trajectories around obstacles using
efficient convex optimization [4, 5, 6, 7, 8]. These motion planners have
demonstrated great potential, and performance frequently superior to widely
used sampling-based methods [5]. However, the decomposition of the environment
that they require is a daunting task, often demanding a substantial degree of
human supervision. This is largely due to the fact that the collision-free
subset of a robot’s configuration space ($\mathcal{C}^{\mathrm{free}}$) is
intractable to describe analytically [9, §3], even if the robot’s task space
is relatively simple (see Figure 1 for an example).

In [10], the authors proposed an algorithm called IRIS to quickly compute
large convex subsets of $\mathcal{C}^{\mathrm{free}}$. This method takes as
input a “seed” configuration and inflates a large polytopic region around this
seed point using convex optimization. While originally limited to the case of
convex obstacles, the IRIS algorithm has been recently extended to handle
nonconvex obstacles and the complicated configuration spaces of multi-link
robot manipulators [11, 12, 13]. The polytopes generated with these algorithms
have been used with great success for motion planning in high dimensions [5,
6, 7]. Nevertheless, seeding collections of these regions so that they cover
diverse areas of $\mathcal{C}^{\mathrm{free}}$ remains a challenge: manual
seeding is tedious, and naively growing regions around randomly chosen
configurations leads to very inefficient decompositions.

Figure 1: The collision-free configuration space of a simple robot is
decomposed into 7 polytopes, achieving around 92% coverage. Left: Robot with
with 3 revolute joints $q_{1}$ to $q_{3}$. Center: Visualization of the full
collision-free configuration space $\mathcal{C}^{\mathrm{free}}$, given by the
interior of the green mesh. Right: Approximate convex cover of
$\mathcal{C}^{\mathrm{free}}$ generated with the proposed method. See also:
Interactive visualization, video.  Figure 2: Sketch of the proposed algorithm
on a simple example. First four figures: Samples are drawn uniformly from
$\mathcal{C}^{\mathrm{free}}$ to build a visibility graph. The visibility
graph is decomposed into five cliques. The principal directions and locations
of the cliques are used to direct a region-inflation algorithm. Remaining two
figures: This process is repeated until sufficient coverage is obtained by
drawing new samples from the remaining free space, and repeating the previous
steps.

In this paper, we propose an efficient method for the approximate
decomposition of robot configuration spaces into few large convex sets,
without any human supervision. A guiding illustration is shown in Figure 2.

Similar to some motion-planning algorithms [14, 15], our method constructs a
visibility graph by sampling points in $\mathcal{C}^{\mathrm{free}}$. The
vertices of this graph are collision-free samples, and the edges connect pairs
of points with mutual line of sight. The visibility graph contains rich
information about the geometry of $\mathcal{C}^{\mathrm{free}}$. In
particular, our key observation is that, as the number of samples grows, fully
connected subgraphs of this visibility graph (so-called “cliques”) tend to
represent better and better approximations of collision-free convex sets in
the underlying configuration space. Our approach is to decompose the
visibility graph into a small collection of large cliques. We then
circumscribe the points in each clique with an ellipsoid by solving a convex-
optimization problem. The center and the principal directions of these
ellipsoids are subsequently used to initialize an inflation algorithm
analogous to IRIS.

Through a large variety of experiments, we show that our algorithm outperforms
the approaches that have been previously used in the literature to seed and
inflate convex regions; both in terms of runtimes and number of regions used
to cover the space. As an example, for a robot arm with seven degrees of
freedom, and many task-space obstacles, our method requires approximately $46$
regions and one hour of computations to cover $70\%$ of
$\mathcal{C}^{\mathrm{free}}$. Whereas the approach outlined in [12] requires
ten times more regions and is ten times slower.

##  II Related Works

Finding the minimum convex cover of a set is a hard problem; even in the case
of two-dimensional polygons, the problem is hard to solve exactly and
approximately.111More formally, the problem is $\exists\mathbb{R}$-complete
[16], and therefore NP-hard, as well as APX-hard [17]. Nonetheless, finding
low-cardinality convex covers of high-dimensional nonconvex spaces (both
polygonal and non-polygonal) remains a problem of practical importance, for
which a variety of approximate algorithms have been devised. In the following,
we group these algorithms into two categories: ones that require explicit
(e.g., analytic) descriptions of $\mathcal{C}^{\mathrm{free}}$, and ones that
only use implicit descriptions of $\mathcal{C}^{\mathrm{free}}$ (and are hence
suitable for most complex configuration spaces). In this literature review, we
particularly focus on IRIS algorithms, due to their efficient scaling to high
dimensions.

###  II-A Algorithms Requiring Explicit Obstacle Descriptions

The recent work [18] constructs low-cardinality convex covers of two-
dimensional polygons by first decomposing them into small convex pieces (for
example triangles), and then by joining the pieces based on a small clique
cover of an approximated set-to-set visibility graph. In three dimensions the
approach in [19] can be used to decompose $\mathcal{C}^{\mathrm{free}}$ into
sets that are approximately convex, provided that a triangle mesh description
is available. Finally, if the configuration-space obstacles are explicitly
described as convex sets, the original IRIS algorithm from [10] can be used to
inflate a large polytope in $\mathcal{C}^{\mathrm{free}}$ around a specified
seed point in arbitrary dimensions.

###  II-B Algorithms Allowing Implicit Obstacle Descriptions

Most commonly, the explicit descriptions of the obstacles are only available
in task space; while the collision-free configuration space
$\mathcal{C}^{\mathrm{free}}$ is defined implicitly through the robot’s
inverse kinematics, and intractable to describe analytically [9, §3].
Obtaining convex decompositions in this setting has been the focus of multiple
works. The method in [20] uses visibility graphs and kernels to compute convex
decompositions of three-dimensional spaces via sample-based collision
checking. However, that method represents the convex sets through their
vertices, and is inefficient in higher dimensions. In the family of IRIS
algorithms, two methods can deal with implicit descriptions of
$\mathcal{C}^{\mathrm{free}}$: IRIS-NP [12] and C-IRIS [11, 13]. The former
extends the original IRIS method [10] to arbitrary configuration spaces using
nonlinear programming and inflates polytopes that are collision-free with high
probability. The latter grows polytopes that are rigorously _certified_ to be
collision-free using a rational reparametrization of the configuration space
and sums-of-squares programming.

##  III Convex Covers, Visibility, and Cliques

In this section, we formally define our main problem: approximating the free
configuration space $\mathcal{C}^{\mathrm{free}}$ with a low-cardinality
collection of convex sets. We also briefly review the main technical tools
that we will use in the development of our algorithm, namely, visibility
graphs and clique covers.

###  III-A Problem Statement

Let $\mathcal{C}^{\mathrm{free}}\subseteq\mathbb{R}^{n}$ be the collision-free
subset of an $n$-dimensional configuration space, which we assume to have a
well-defined finite volume. Let also $\alpha$ be a constant in the interval
$(0,1]$.

###### Definition 1

An _$\alpha$ -approximate convex cover_ of $\mathcal{C}^{\mathrm{free}}$ is a
collection of potentially overlapping convex sets
$\mathcal{R}_{1},\ldots,\mathcal{R}_{N}\subseteq\mathcal{C}^{\mathrm{free}}$
whose union covers at least an $\alpha$-fraction of the volume of
$\mathcal{C}^{\mathrm{free}}$:

$\mathop{\bf
vol}\left(\bigcup_{i=1}^{N}\mathcal{R}_{i}\right)\geq\alpha\mathop{\bf
vol}\left(\mathcal{C}^{\mathrm{free}}\right).$

Our problem is to find an $\alpha$-approximate convex cover of minimum
cardinality $N$.

Problem:  Min $\alpha$-ApproxConvexCover
$\displaystyle\mathop{\text{minimize}}N\mathop{\text{subject to}}$
$\displaystyle\mathop{\bf
vol}\left(\bigcup_{i=1}^{N}\mathcal{R}_{i}\right)\geq\alpha\mathop{\bf
vol}\left(\mathcal{C}^{\mathrm{free}}\right),$
$\displaystyle\mathcal{R}_{i}\subseteq\mathcal{C}^{\mathrm{free}},\quad\forall
i=1,\dots,N.$

In practice, we are interested in solving this problem for values of $\alpha$
that are sufficiently high to accomplish a task of interest, such as
collision-free motion planning, but not so large that the cardinality $N$
grows unreasonably. Indeed, when $\alpha=1$, Min $\alpha$-ApproxConvexCover
might not even have a finite solution.222For $\alpha<1$ a finite solution is
guaranteed to exist. This can be seen by approximating the volume of
$\mathcal{C}^{\mathrm{free}}$ over finite hyperrectangle partitions with a
lower Darboux integral.

###  III-B Visibility Graphs and Clique Covers

Our algorithm is based on the idea that clusters of points that see each other
can approximate convex subsets of $\mathcal{C}^{\mathrm{free}}$. Here, we
formally define the notion of visibility as well as introduce some formal
tools from graph theory to guide the development of our algorithm.

We begin by defining visibility in $\mathcal{C}^{\mathrm{free}}$.

###### Definition 2

Two points $q,q^{\prime}\in\mathcal{C}^{\mathrm{free}}$ are said to _see each
other_ if the entire line connecting them is collision-free:
$tq+(1-t)q^{\prime}\in\mathcal{C}^{\mathrm{free}}$ for all $t\in[0,1]$. Notice
that this definition is symmetric in $q$ and $q^{\prime}$.

We are now ready to define the visibility graph of a set of collision-free
configurations.

###### Definition 3

The _visibility graph_ of a set of points
$q_{1},\dots,q_{K}\in\mathcal{C}^{\mathrm{free}}$ is an undirected graph
$\mathcal{G}=(\mathcal{V},\mathcal{E})$ with vertices
$\mathcal{V}=\\{1,\ldots,K\\}$, and with an edge $\\{i,j\\}\in\mathcal{E}$ for
every pair of distinct points $q_{i}$ and $q_{j}$ that see each other.

We show an example of a visibility graph in Figure 2. We note that clusters of
points that can all see each other form a clique in the visibility graph.

###### Definition 4

Let $\mathcal{G}=(\mathcal{V},\mathcal{E})$ be an undirected graph. A _clique_
$\mathcal{K}$ is a subset of $\mathcal{V}$ where every pair of vertices is
connected by an edge.

Note that if a cluster of configurations can be placed in the same convex set,
then these configurations must form a clique in the visibility graph. The
second panel in Figure 2 highlights a collection of five cliques in the
visibility graph that have this property. These five cliques form what is
called a clique cover; which resembles a discrete analog of a convex cover.

###### Definition 5

A collection $\mathcal{T}$ of cliques $\mathcal{K}_{1},\dots,\mathcal{K}_{N}$
is a _clique cover_ of a graph $\mathcal{G}$ if every vertex in the graph is
contained in at least one clique.

A natural discrete counterpart of the minimum convex cover is the
MinCliqueCover problem, where we look for the minimum number of cliques $N$
required to cover a graph. Our observation is that, as the number of samples
in the visibility graph increases, a minimum clique cover typically does an
increasingly better job of approximating a minimum convex cover. Limitations
of this analogy are discussed in §VI.

MinCliqueCover is NP-complete [21]. There exist heuristics, such as [22], that
attempt to solve MinCliqueCover directly. Alternatively, one can greedily
construct a clique cover by repeatedly eliminating the largest clique. The
problem of finding the largest clique in a graph is called MaxClique. Even
though this problem is also NP-complete [21], it is often substantially faster
to solve in practice. We found the latter approach with exact solutions of
MaxClique to perform particularly well on our problem instances.

##  IV Algorithm

We now present our Visibility Clique Cover (VCC) algorithm, which consists of
four main steps. First, we randomly sample a collection of points in
$\mathcal{C}^{\mathrm{free}}$ and construct their visibility graph. Second, we
compute an approximate clique cover of the graph. Third, we summarize the
geometric information of each clique using an ellipsoid. Fourth, we use these
ellipsoids to initialize a polytope-inflation algorithm analogous to IRIS.
This process is repeated until the generated set of polytopes $\mathcal{R}$
covers a given fraction $\alpha$ of $\mathcal{C}^{\mathrm{free}}$. This
procedure is summarized in Algorithm 1, and the remainder of this section
details the individual steps.

Input :

$\alpha$: coverage threshold

$K$: number of samples per iteration

$s_{\text{min}}$: minimum clique size

Output :

$\mathcal{R}$: set of convex polytopes approximating
$\mathcal{C}^{\mathrm{free}}$

Algorithm: __

$\mathcal{R}\leftarrow\emptyset$

while  _$\textsc{CheckCoverage}(\mathcal{R})\leq\alpha$_ do

$\mathcal{G}~{}\leftarrow\textsc{SampleVisibilityGraph}(K,\mathcal{R})$

$\mathcal{T}\leftarrow\textsc{TruncatedCliqueCover}(\mathcal{G},s_{\text{min}})$

$\mathcal{B}\leftarrow\textsc{MinVolumeEllipsoids}(\mathcal{T})$

$\mathcal{R}\leftarrow\mathcal{R}\cup\textsc{InflatePolytopes}(\mathcal{B})$

end while

return $\mathcal{R}$

Algorithm 1 VisibilityCliqueInflation

###  IV-A Sampling the Visibility Graph

At the beginning of every iteration of VCC, the subroutine
SampleVisibilityGraph samples $K$ configurations uniformly at random from the
portion of $\mathcal{C}^{\mathrm{free}}$ that is not already covered by the
polytopes in $\mathcal{R}$. Then, it constructs the visibility graph
$\mathcal{G}=(\mathcal{V},\mathcal{E})$, by checking for collisions along the
line segments connecting each pair of sampled configurations. Currently, this
is performed using sampling-based collision checkers. Exact visibility
checking is possible using methods such as [23], [24, §5.3.4], or [13].

###  IV-B Truncated Clique Cover

In the subroutine TruncatedCliqueCover we approximately cover the visibility
graph with a collection of cliques, each of which contain at least $s_{\min}$
vertices. We construct this approximate cover $\mathcal{T}$ greedily, by
solving a sequence of MaxClique problems. Each instance of MaxClique is
formulated as an integer linear program

$\displaystyle\mathop{\text{maximize}}\quad$
$\displaystyle\sum_{i=1}^{K}b_{i}$ (1a) $\displaystyle\mathop{\text{subject
to}}\quad$ $\displaystyle b_{i}+b_{j}\leq 1,$
$\displaystyle\forall\\{i,j\\}\in\bar{}\mathcal{E},$ (1b) $\displaystyle
b_{i}\in\\{0,1\\},$ $\displaystyle\forall i=1,\ldots,K.$ (1c)

A binary decision variable $b_{i}$ is added for each vertex. The role of this
variable is to take unit value if and only if vertex $i$ is included in the
clique. The set $\bar{}\mathcal{E}$ contains all the pairs of vertices
$\\{i,j\\}$ such that $i\neq j$ and $\\{i,j\\}\notin\mathcal{E}$. Therefore
the first constraint ensures that two vertices are selected only if they share
an edge.

After solving the integer program (1), the clique found is removed from the
graph and added to the clique cover. Since small cliques are not informative,
we stop this iterative process when the largest clique left in the graph is
smaller than a given threshold $s_{\min}$. For this reason, our clique covers
will be truncated, i.e., generally, not all vertices will be contained in one
of the cliques.

###  IV-C Summarizing Cliques with Ellipsoids

In the subroutine MinVolumeEllipsoids, we solve a semidefinite program to
enclose each clique with an ellipsoid of minimum volume [25, §8.4.1]. This
collection $\mathcal{B}$ of ellipsoids allows us to summarize the geometry of
each clique with a point and a set of principal directions, which are then
used to initialize the region-inflation algorithm. For the upcoming
computations, it is necessary that the center of each ellipsoid is not in
collision; if this is not the case, we recenter the ellipsoid around the
vertex in the clique that is closest to its center.

###  IV-D Inflating Polytopes

Figure 3: The growth direction of an IRIS region can be guided by the initial
distance metric. IRIS is initialized with three ellipsoids with same center
but different principal axes, resulting in polytopes that cover different
portions of $\mathcal{C}^{\mathrm{free}}$.

In the last step of VCC, the subroutine InflatePolytopes inflates a large
collision-free polytope around the center of each ellipsoid induced by a
clique.

Let us initially assume that the obstacles are convex. Consider a single
ellipsoid. Using convex optimization we compute the point in each obstacle
that is closest to the center of the ellipsoid, according to the distance
metric induced by the ellipsoid. These points anchor separating hyperplanes
between the ellipsoid center and the obstacles, which form a polytope of
obstacle-free space. These steps are repeated for each ellipsoid (i.e., for
each clique) to obtain a collection of collision-free polytopes that we add to
the set $\mathcal{R}$.

(a) Village

(b) 3DOF Flipper

(c) 5DOF UR5

(d) 7DOF IIWA

Domain | Village | 3DOF Flipper | 5DOF UR5 | 7DOF IIWA
---|---|---|---|---
Algorithm | IOS | VCC (ours) | IOS | VCC (ours) | IOS | VCC (ours) | IOS | VCC (ours)
# regions $|\mathcal{R}|$ | 198.0$\pm$13.6 |  93.9$\pm$7.5 | 10.4$\pm$1.9 |  6.7$\pm$0.5 | 90.5$\pm$18.8 |  35.1$\pm$1.6 | 482.6$\pm$83.8 |  46.3$\pm$4.5
runtime [s] | 2.7e3$\pm$2.9e2 |  1.7e3$\pm$3.4e2 | 2.0e2$\pm$2.9e1 |  4.9e1$\pm$7.1 | 1.12e4$\pm$2.1e3 |  5.5e2$\pm$6.6e1 | 8.9e4$\pm$1.9e4 |  5.7e3$\pm$ 1.7e3
coverage threshold $\alpha$ | 0.8 | 0.8 | 0.9 | 0.9 | 0.75 | 0.75 | 0.7 | 0.7
# visibility vertices $K$ | - | 500 | - | 500 | - | 1000 | - | 1500
min. clique size $s_{\text{min}}$ | - | 10 | - | 10 | - | 10 | - | 20

TABLE I: Comparison of our Visibility Clique Cover (VCC) algorithm to the
Iterative Obstacle Seeding (IOS) method from [12], across four different
environments. All experiments are repeated ten times. The numbers in the first
two rows indicate the mean and the empirical standard deviation over the
trials. We observe that VCC achieves the given coverage targets with 1.6 to
10.4 times fewer regions, and between 1.4 to 20 times faster than IOS. The
environment names are linked to interactive visualizations.

These computations correspond to a single iteration of the IRIS algorithm
[10], and ensure that the largest uniformly-scaled, collision-free version of
the ellipsoid is contained in the resulting polytope. Figure 3 illustrates how
the initial metric is fundamental in guiding the shape of the regions
generated by IRIS. Traditionally, the IRIS algorithm was initialized with an
uninformed ellipsoidal metric and needed to run for multiple (expensive)
iterations in order to expand and cover a larger volume of space; in VCC we
require only a single IRIS iteration.

When the obstacles are not convex, we run one iteration of the nonlinear
programming variant IRIS-NP [12] instead. Alternatively, C-IRIS [11, 13] could
be employed to obtain certifiably collision-free regions.

###  IV-E Convergence Check

The subroutine CheckCoverage estimates the fraction of
$\mathcal{C}^{\mathrm{free}}$ covered by the regions in $\mathcal{R}$, and
terminates our algorithm if this value exceeds the threshold $\alpha$.
Computing this fraction exactly is impractical, and so we resort to randomized
methods. The coverage is estimated by drawing a large number of $M$ samples in
$\mathcal{C}^{\mathrm{free}}$, and computing the ratio of samples that land in
at least one of the regions in $\mathcal{R}$. More sophisticated checks, such
as one-sided Bernoulli hypothesis testing, are possible.

###  IV-F Completeness

Analogous to [26, 27], VCC is probabilistically complete under mild
assumptions. This means as the number of iterations goes to infinity, the
probability of completely covering $\mathcal{C}^{\mathrm{free}}$ goes to one.

##  V Experiments

As there are no direct baseline methods, we compare our VCC algorithm against
an extension of the method in [12, §III.D]. The natural extension of this
approach is to iteratively grow polytopes around uniformly sampled points from
the uncovered free space using IRIS, while treating previously computed
regions as obstacles. This process is repeated until the desired coverage is
met. We call this approach Iterative Obstacle Seeding (IOS).

In the following, all experiments are implemented in Drake [28], and all
computations are performed on a single desktop computer with an Intel Core
i9-10850K CPU and 32 Gb of RAM. We solve all MaxClique instances to global
optimality using Gurobi [29].

We evaluate VCC and IOS on four environments: Village, 3DOF Flipper, 5DOF UR5,
and 7DOF IIWA. The dimension $n$ of these environments ranges from $3$ to $7$.
For each environment, we run the two algorithms ten times and report their
performance in Table I. The Village environment from [8] contains only convex
obstacles and is compatible with the original IRIS [10]. All other examples
involve the configuration spaces of robotic manipulators and therefore IRIS-NP
[12] is employed. The number of samples used in the convergence check in
Algorithm 1 is $M=5000$.

VCC meets the required coverage threshold with significantly fewer regions and
in a substantially shorter amount of time. Notably, in the most challenging
benchmark, 7DOF IIWA, VCC requires 10 times fewer regions and meets the
required coverage of $70\%$ around 16 times faster. This substantial speedup
can be attributed to two key factors. First, the region inflation in VCC is
parallelized. Second, IRIS is the most computationally expensive step. While
VCC only requires a single iteration of IRIS per region, IOS can require up to
around five IRIS iterations per region, due to the initialization with a
potentially uninformative spherical metric.

##  VI Limitations of Approximating Convex Sets with Cliques

Despite the strong performance of our algorithm, the hardness of solving Min
$\alpha$-ApproxConvexCover makes it possible to construct simple examples that
highlight pitfalls of our heuristic approach. In this section, we discuss
holes in $\mathcal{C}^{\mathrm{free}}$ which leads to one such pitfall that is
particularly insightful.

While every convex set in $\mathcal{C}^{\mathrm{free}}$ naturally corresponds
to a clique in the visibility graph, a clique in the visibility graph does not
necessarily correspond to a convex set in $\mathcal{C}^{\mathrm{free}}$. The
convex hull of a clique can enclose holes. This problem persists even if we
sample arbitrarily dense visibility graphs, and if we restrict the analysis to
maximum cliques. A visual proof is shown in Figure 5, which illustrates a
triangular configuration space with a triangular hole of size $\varepsilon$.
As $\varepsilon$ goes to zero, the largest convex subset of
$\mathcal{C}^{\mathrm{free}}$ is the green trapezoid, whose area approaches
$5/9$ of the total area of $\mathcal{C}^{\mathrm{free}}$. On the other hand, a
larger subset of mutually visible points is given by the union of the three
red parallelograms, whose area approaches $6/9$ of the total area. Therefore,
if configurations are sampled uniformly at random, an optimum solution of
MaxClique will almost surely enclose a hole as the number of samples goes to
infinity. A similar construction can be used to show an analogous discrepancy
between the minimum convex cover of $\mathcal{C}^{\mathrm{free}}$ and
MinCliqueCover.

In principle, this problem can be addressed by solving a modified version of
MaxClique that better captures the notion of a convex set. In short, we
require that every vertex of the visibility graph that is contained in the
convex hull of the maximum clique must be a member of the clique. For an
infinitely dense visibility graph, this ensures that the maximum clique cannot
enclose holes.

$\varepsilon$1 Figure 5: Maximum cliques of infinitely dense visibility graphs
can enclose holes, and do not necessarily correspond to collision-free convex
sets. The largest collision-free convex region (green trapezoid) has a smaller
area than the union of red parallelograms when $0<\varepsilon\leq
1-\sqrt{5/6}$. In this case, the convex hull of the maximum clique must
enclose the hole.

We enforce the contrapositive of the latter condition through linear
constraints that separate all non-clique members $q_{i}$, with
$i\in\\{1,\ldots,K\\}$, from the points in the clique with a hyperplane
$\mathcal{H}_{i}=\\{q\mid c_{i}^{T}q+d_{i}=0\\}$, parameterized by the
decision variables $(c_{i},d_{i})\in\mathbb{R}^{n+1}$.

The resulting mixed-integer linear optimization extends (1) by adding the
additional separation constraints

$\displaystyle c_{i}^{T}q_{i}+d_{i}$ $\displaystyle\geq 1-b_{i},\quad$
$\displaystyle\forall~{}i=1,\dots,K,$ (2a) $\displaystyle
c_{i}^{T}q_{j}+d_{i}$ $\displaystyle\leq M(1-b_{j}),\quad$
$\displaystyle\forall i,j=1,\dots,K,$ (2b)

where $M$ is a large enough constant (e.g. the maximum distance between all
pairs of vertices).

When the point $q_{i}$ is not in the clique ($b_{i}=0$) and the point $q_{j}$
is in the clique ($b_{j}=1$), these constraints read $c_{i}^{T}q_{i}+d_{i}\geq
1$ and $c_{i}^{T}q_{j}+d_{i}\leq 0$. Therefore, the hyperplane
$\mathcal{H}_{i}$ separates $q_{i}$ from $q_{j}$. On the other hand, these
constraints are seen to be redundant for all other possible values of the
binaries $b_{i}$ and $b_{j}$. In Figure 6, we demonstrate how this extension
prevents a maximum clique from enclosing holes in a finite-sample regime.

In practice, solving (1) with constraints (2) is too expensive for the
problems in Table I. Nonetheless, we observed that in the first MaxClique
problem (1) only an average of $0.1\%$ of the vertices were excluded from the
clique that could not be separated from it with a hyperplane. Subsequent
cliques, and improvements to the formulation (2), demand a more nuanced
discussion and are subject to future work.

##  VII Conclusions

We have proposed an algorithm for approximately decomposing complex
configuration spaces into small collections of polytopes. Our algorithm uses
clique covers of visibility graphs as an effective heuristic for obtaining
local information about $\mathcal{C}^{\mathrm{free}}$, and for seeding a
region-inflation algorithm. The parallels between convex sets in
$\mathcal{C}^{\mathrm{free}}$ and cliques in visibility graphs have also been
discussed. Our experiments demonstrate that VCC reliably finds approximate
convex covers of $\mathcal{C}^{\mathrm{free}}$ with fewer regions and in less
time than previous approaches.

Figure 6: A visibility graph with 100 random vertices in the triangular
environment from Figure 5. Solving the maximum clique problem (1) with the
additional constraints (2) yields a clique with 56 vertices (shown in green),
that closely approximates the corresponding convex set in Figure 5. Solving
only problem (1), yields a clique with 63 vertices (shown in red) that,
however, encloses the central hole.

## Acknowledgements

The authors would like to thank Steven M. LaValle, Thomas Cohn, Annan Zhang,
and Fabian Dickhardt for the many fruitful discussions.

## References

  * [1] C. M. Bishop and N. M. Nasrabadi, _Pattern recognition and machine learning_. Springer, 2006, vol. 4, no. 4.
  * [2] M. Pharr, C. Kolb, R. Gershbein, and P. Hanrahan, “Rendering complex scenes with memory-coherent ray tracing,” in _Proceedings of the 24th annual conference on Computer graphics and interactive techniques_ , 1997, pp. 101–108.
  * [3] E. G. Gilbert, D. W. Johnson, and S. S. Keerthi, “A fast procedure for computing the distance between complex objects in three-dimensional space,” _IEEE Journal on Robotics and Automation_ , vol. 4, no. 2, pp. 193–203, 1988\.
  * [4] T. Marcucci, J. Umenberger, P. A. Parrilo, and R. Tedrake, “Shortest paths in graphs of convex sets,” _arXiv preprint arXiv:2101.11565_ , 2021.
  * [5] T. Marcucci, M. Petersen, D. von Wrangel, and R. Tedrake, “Motion planning around obstacles with convex optimization,” _arXiv preprint arXiv:2205.04422_ , 2022.
  * [6] T. Cohn, M. Petersen, M. Simchowitz, and R. Tedrake, “Non-euclidean motion planning with graphs of geodesically-convex sets,” _arXiv preprint arXiv:2305.06341_ , 2023.
  * [7] V. Kurtz and H. Lin, “Temporal logic motion planning with convex optimization via graphs of convex sets,” _arXiv preprint arXiv:2301.07773_ , 2023.
  * [8] T. Marcucci, P. Nobel, R. Tedrake, and S. Boyd, “Fast path planning through large collections of safe boxes,” _arXiv preprint arXiv:2305.01072_ , 2023\.
  * [9] J.-C. Latombe, _Robot motion planning_. Springer Science & Business Media, 2012, vol. 124.
  * [10] R. Deits and R. Tedrake, “Computing large convex regions of obstacle-free space through semidefinite programming,” in _Algorithmic Foundations of Robotics XI_. Springer, 2015, pp. 109–124.
  * [11] A. Amice, H. Dai, P. Werner, A. Zhang, and R. Tedrake, “Finding and optimizing certified, collision-free regions in configuration space for robot manipulators,” in _Algorithmic Foundations of Robotics XV: Proceedings of the Fifteenth Workshop on the Algorithmic Foundations of Robotics_. Springer, 2022, pp. 328–348.
  * [12] M. Petersen and R. Tedrake, “Growing convex collision-free regions in configuration space using nonlinear programming,” _arXiv preprint arXiv:2303.14737_ , 2023.
  * [13] H. Dai, A. Amice, P. Werner, A. Zhang, and R. Tedrake, “Certified polyhedral decompositions of collision-free configuration space,” _arXiv preprint arXiv:2302.12219_ , 2023.
  * [14] T. Siméon, J.-P. Laumond, and C. Nissoux, “Visibility-based probabilistic roadmaps for motion planning,” _Advanced Robotics_ , vol. 14, no. 6, pp. 477–493, 2000.
  * [15] T. Lozano-Pérez and M. A. Wesley, “An algorithm for planning collision-free paths among polyhedral obstacles,” _Communications of the ACM_ , vol. 22, no. 10, pp. 560–570, 1979.
  * [16] M. Abrahamsen, “Covering polygons is even harder,” in _2021 IEEE 62nd Annual Symposium on Foundations of Computer Science (FOCS)_. IEEE, 2022, pp. 375–386.
  * [17] S. J. Eidenbenz and P. Widmayer, “An approximation algorithm for minimum convex cover with logarithmic performance guarantee,” _SIAM Journal on Computing_ , vol. 32, no. 3, 2003.
  * [18] M. Abrahamsen, W. Bille Meyling, and A. Nusser, “Constructing concise convex covers via clique covers (cg challenge),” in _39th International Symposium on Computational Geometry (SoCG 2023)_. Schloss Dagstuhl-Leibniz-Zentrum für Informatik, 2023.
  * [19] K. Mamou and F. Ghorbel, “A simple and efficient approach for 3d mesh approximate convex decomposition,” in _2009 16th IEEE international conference on image processing (ICIP)_. IEEE, 2009, pp. 3501–3504.
  * [20] A. Sarmientoy, R. Murrieta-Cidz, and S. Hutchinsony, “A sample-based convex cover for rapidly finding an object in a 3-d environment,” in _Proceedings of the 2005 IEEE International Conference on Robotics and Automation_. IEEE, 2005, pp. 3486–3491.
  * [21] R. M. Karp, _Reducibility among combinatorial problems_. Springer, 2010.
  * [22] D. Strash and L. Thompson, “Effective data reduction for the vertex clique cover problem,” in _2022 Proceedings of the Symposium on Algorithm Engineering and Experiments (ALENEX)_. SIAM, 2022, pp. 41–53.
  * [23] F. Schwarzer, M. Saha, and J.-C. Latombe, “Exact collision checking of robot paths,” _Algorithmic foundations of robotics V_ , pp. 25–41, 2004.
  * [24] S. M. LaValle, _Planning algorithms_. Cambridge university press, 2006.
  * [25] S. P. Boyd and L. Vandenberghe, _Convex optimization_. Cambridge university press, 2004.
  * [26] J. J. Kuffner and S. M. LaValle, “RRT-connect: An efficient approach to single-query path planning,” in _Proceedings 2000 ICRA. Millennium Conference. IEEE International Conference on Robotics and Automation. Symposia Proceedings (Cat. No. 00CH37065)_ , vol. 2. IEEE, 2000, pp. 995–1001.
  * [27] L. E. Kavraki, P. Svestka, J.-C. Latombe, and M. H. Overmars, “Probabilistic roadmaps for path planning in high-dimensional configuration spaces,” _IEEE transactions on Robotics and Automation_ , vol. 12, no. 4, pp. 566–580, 1996.
  * [28] R. Tedrake and the Drake Development Team, “Drake: Model-based design and verification for robotics,” 2019. [Online]. Available: https://drake.mit.edu
  * [29] Gurobi Optimization, LLC, “Gurobi Optimizer Reference Manual,” 2023.


# AMICI: High-Performance Sensitivity Analysis for Large Ordinary Differential
Equation Models

Fabian Fröhlich 1,∗, Daniel Weindl 2, Yannik Schälte 2,3,
Dilan Pathirana 4, Łukasz Paszkowski 5, Glenn Terje Lines 5,
Paul Stapor 2,3 and Jan Hasenauer 2,3,4,∗

(December 2020)

###### Abstract

Summary: Ordinary differential equation models facilitate the understanding of
cellular signal transduction and other biological processes. However, for
large and comprehensive models, the computational cost of simulating or
calibrating can be limiting. AMICI is a modular toolbox implemented in
C++/Python/MATLAB that provides efficient simulation and sensitivity analysis
routines tailored for scalable, gradient-based parameter estimation and
uncertainty quantification.
Availability: AMICI is published under the permissive BSD-3-Clause license
with source code publicly available on https://github.com/AMICI-dev/AMICI.
Citeable releases are archived on Zenodo.
Contact<EMAIL_ADDRESS><EMAIL_ADDRESS>
Supplementary information: Supplementary information is available at arXiv
online.

1Department of Systems Biology, Harvard Medical School, Boston, MA 02115, USA,
2Institute of Computational Biology, Helmholtz Zentrum München – German
Research Center for Environmental Health, 85764 Neuherberg, Germany,
3Center for Mathematics, Technische Universität München, 85748 Garching,
Germany,
4Faculty of Mathematics and Natural Sciences, University of Bonn, 53113 Bonn,
Germany, and
5Simula Research, 1325 Lysaker, Norway.
∗To whom correspondence should be addressed.

##  1 Introduction

Ordinary Differential Equation (ODE) models are widely used in systems
biology, for example, to elucidate dynamic processes and to predict response
to perturbations. Model parameters have to be inferred from data, which can be
computationally intensive as thousands of simulations may be required. This is
challenging for large models, with many state variables and parameters, where
simulations take seconds to minutes (Fröhlich et al.,, 2018). Inference can
benefit from accurate sensitivities (Villaverde et al.,, 2018), i.e.,
derivatives of model outputs with respect to model parameters. Accurate
sensitivities can be computed using forward, adjoint or steady-state
approaches, but benefit from symbolic derivatives of model equations, which
are labor-intensive and error-prone to compute manually.

The recent surge of genome-scale perturbation data as well as respective
comprehensive models has increased the demand for methods for scalable
sensitivity computation for ODE models. To address this demand, we introduce
AMICI, a high-performance simulation and sensitivity analysis library. AMICI
is implemented in C++ and Python, and provides model import from widely used
formats such as the Systems Biology Markup Language (SBML) (Hucka et al.,,
2003), BioNetGen Language (BNGL) (Harris et al.,, 2016) and Kappa (Boutillier
et al.,, 2018), and generates high-performance-computing (HPC) ready modules.
These modules provide model-specific simulation and sensitivity computation
routines, which can be accessed from Python, C++, and MATLAB. For parameter
estimation problems specified in PEtab (Schmiester et al.,, 2020), AMICI can
evaluate the objective function and its gradient.

##  2 Methods

For high simulation performance, AMICI reads models from high-level formats,
then translates the model and symbolically derived expressions to C++ code. As
symbolic processing can be computationally intensive, AMICI symbolically only
computes partial derivatives, total derivatives are computed through (sparse)
matrix multiplication and addition at runtime. This expedites model
compilation and simulation.

To simulate stiff, large models, efficient linear solvers are crucial. AMICI
features several direct dense, direct sparse and implicit solvers. For
sensitivity analysis, AMICI implements forward, adjoint (Fröhlich et al.,,
2017), and steady-state (Lines et al.,, 2019) methods and combinations thereof
(see Supplementary Information).

##  3 Implementation

The AMICI library is implemented in C++14 and relies on SUNDIALS (Hindmarsh et
al.,, 2005) and SuiteSparse (Davis and Palamadai Natarajan,, 2010) for
simulation and sensitivity computation. For matrix multiplication, a CBLAS-
compatible library is used (e.g. ATLAS, OpenBLAS, Intel MKL).

Model import is implemented in Python and supports several widely used formats
(Figure 1). SBML import is implemented using libSBML (Bornstein et al.,,
2008), Kappa and BNGL import is implemented via PySB (Lopez et al.,, 2013).
Symbolic processing is performed using SymPy (Meurer et al.,, 2017).
Compilation has been tested with GCC, Clang, Intel and MinGW on Linux, MacOS
and Windows platforms. For interoperability, AMICI provides CMake files and
Dockerfiles.

AMICI ships with a MATLAB interface and implements a Python interface using
SWIG, which can be extended to other languages. Simulation and compilation are
highly configurable with the API documented on Read the Docs and in example
notebooks. For HPC-readiness, AMICI implements OpenMP parallelization over
experimental conditions and is thread-safe.

Simulation and sensitivity analysis rely on intricate theory that is complex,
error-prone to implement. Accordingly, we implemented an extensive continuous
integration test pipeline. Simulation results are verified for the SBML
semantic test suite, where AMICI passes 862 out of 1780 tests (appropriate
error messages to indicate unsupported features for remaining tests).
Simulation results are compared against PySB simulations for 17 BNGL
validation models and examples, where AMICI passes all tests. Sensitivity
results are checked by regression tests covering, e.g. forward, adjoint and
steady-state sensitivities. Performance tests check computation time for model
import, simulation and sensitivity analysis of a large model. Unit and
documentation tests, static code analysis, and memory leak checks are
included, and code review is enforced for all contributions.

Figure 1: Outline of model compilation in AMICI. Models, specified in SBML or
PySB, are imported in Python. This generates efficient C++ code for model
equations and simulation bindings for Python, C++ and MATLAB.

##  4 Discussion

There is a rich ecosystem of tools for model simulation. To avoid duplication,
we ensured good interoperability with other tools: AMICI does not provide a
model development environment, but permits model import from standard formats.
Similarly, AMICI is not part of an integrated parameter estimation framework,
but features a flexible, well-documented API. Currently, several parameter
estimation tools can interface AMICI, including pyPESTO (Schälte et al.,,
2020) and parPE (Schmiester et al.,, 2019). This modular design aims at
researchers developing new methods or tools for parameter estimation that
would benefit from state-of-the-art simulation and sensitivity computation
routines. Specifically, adjoint and steady-state sensitivity analysis as well
as sparse linear solvers are currently only supported in a small, disparate
set of tools, which highlights the unique capabilities of AMICI.

AMICI has been in development since 2015 and has so far been used in at least
50 publications, and is continuously being developed by 4 core contributors at
3 different institutions. In the future, we plan to improve SBML support and
extend interoperability with other tools.

## Funding

This work was supported by the European Union’s Horizon 2020 research and
innovation program (CanPathPro; Grant no. 686282; J.H., D.W., P.S., Ł.P.,
G.T.L., F.F.), the Federal Ministry of Education and Research of Germany
(Grant no. 01ZX1916A; D.W. & 01ZX1705A; J.H., & Grant. no. 031L0159C; J.H.),
the German Research Foundation (Grant no. HA7376/1-1; Y.S., Germany’s
Excellence Strategy – EXC-2047/1 – 390685813; D.P.), the Human Frontier
Science Program (Grant no. LT000259/2019-L1; F.F.), the National Institute of
Health (Grant no. U54-CA225088; F.F.), and the Federal Ministry of Economic
Affairs and Energy (Grant no. 16KN074236; D.P.).

## References

  * Bornstein et al., (2008) Bornstein, B. J., Keating, S. M., Jouraku, A., and Hucka, M. (2008).  LibSBML: An API library for SBML.  Bioinformatics, 24(6):880–881.
  * Boutillier et al., (2018) Boutillier, P., Maasha, M., Li, X., Medina-Abarca, H. F., Krivine, J., Feret, J., Cristescu, I., Forbes, A. G., and Fontana, W. (2018).  The Kappa platform for rule-based modeling.  Bioinformatics, 34(13):i583–i592.
  * Davis and Palamadai Natarajan, (2010) Davis, T. A. and Palamadai Natarajan, E. (2010).  Algorithm 907: KLU, a direct sparse solver for circuit simulation problems.  ACM T. Math. Software., 37(3):36.
  * Fröhlich et al., (2018) Fröhlich, F. et al. (2018).  Efficient parameter estimation enables the prediction of drug response using a mechanistic pan-cancer pathway model.  Cell Syst., 7(6):567–579.e6.
  * Fröhlich et al., (2017) Fröhlich, F., Kaltenbacher, B., Theis, F. J., and Hasenauer, J. (2017).  Scalable parameter estimation for genome-scale biochemical reaction networks.  PLoS Comput Biol, 13(1):e1005331.
  * Harris et al., (2016) Harris, L. A., Hogg, J. S., Tapia, J.-J., Sekar, J. A. P., Gupta, S., Korsunsky, I., Arora, A., Barua, D., Sheehan, R. P., and Faeder, J. R. (2016).  BioNetGen 2.2: advances in rule-based modeling.  Bioinformatics, 32(21):3366–3368.
  * Hindmarsh et al., (2005) Hindmarsh, A. C., Brown, P. N., Grant, K. E., Lee, S. L., Serban, R., Shumaker, D. E., and Woodward, C. S. (2005).  SUNDIALS: Suite of Nonlinear and Differential/Algebraic Equation Solvers.  ACM T. Math. Software., 31(3):363–396.
  * Hucka et al., (2003) Hucka, M. et al. (2003).  The systems biology markup language (SBML): A medium for representation and exchange of biochemical network models.  Bioinformatics, 19(4):524–531.
  * Lines et al., (2019) Lines, G. T., Paszkowski, Ł., Schmiester, L., Weindl, D., Stapor, P., and Hasenauer, J. (2019).  Efficient computation of steady states in large-scale ode models of biochemical reaction networks.  IFAC-PapersOnLine, 52(26):32 – 37.  8th Conference on Foundations of Systems Biology in Engineering FOSBE 2019\.
  * Lopez et al., (2013) Lopez, C. F., Muhlich, J. L., Bachman, J. A., and Sorger, P. K. (2013).  Programming biological models in python using pysb.  Mol. Syst. Biol., 9:646.
  * Meurer et al., (2017) Meurer, A. et al. (2017).  Sympy: symbolic computing in python.  PeerJ Computer Science, 3:e103.
  * Schälte et al., (2020) Schälte, Y., Fröhlich, F., Stapor, P., Wang, D., Weindl, D., Schmiester, L., and Hasenauer, J. (2020).  ICB-DCM/pyPESTO: pyPESTO 0.0.11.  https://doi.org/10.5281/zenodo.3715448.
  * Schmiester et al., (2020) Schmiester, L. et al. (2020).  Petab – interoperable specification of parameter estimation problems in systems biology.  arXiv.
  * Schmiester et al., (2019) Schmiester, L., Schälte, Y., Fröhlich, F., Hasenauer, J., and Weindl, D. (2019).  Efficient parameterization of large-scale dynamic models based on relative measurements.  Bioinformatics, 36(2):594–602.
  * Villaverde et al., (2018) Villaverde, A. F., Froehlich, F., Weindl, D., Hasenauer, J., and Banga, J. R. (2018).  Benchmarking optimization methods for parameter estimation in large kinetic models.  Bioinformatics, page bty736.


# A short introduction on Angular momentum of Kerr Blackhole

Sumit Panganti and Siba Prasad Das  Department of Physics, Shivaji University,
Kolhapur-416004, Maharashtra, India
<EMAIL_ADDRESS>
(22 May 2024; 22 May 2024)

###### Abstract

General relativity (GR) predicts the existence of black hole (BH). The
rotating BH called as a Kerr Black hole and GR implies that there is an upper
limit on the angular momentum per mass squared of black holes $\leq 1$, above
which the event horizon of the Kerr BH is not exist. We find the radial
equation for equatorial motion for Kerr BH in terms of the effective
potential. We have shown the effective potential profile for different
rotation parameter ($a$). We find the solution of the radial equation of the
Kerr metric and found the expression of the angular momentum per unit mass
squared, $\tilde{a}=\frac{a}{M}$. We showed the profile of $\tilde{a}$ as a
function of $\frac{r}{M}$. The solution also leads the energy per unit rest
mass ($e$) and we showed its behavior as a function of $\frac{r}{M}$. We
enumerated the maximum values of radius of innermost stable circular orbit
($r_{ISCO}$) for $\tilde{a}=1$.

###### keywords:

Kerr Blackholes; Solutions for Radial Equation; Innermost stable circular
orbit

PACS numbers: 04.70.Bw

##  1 Introduction

The horizon-scale images of black holes (BHs) and their shadows have opened an
unprecedented window onto tests of gravity [1]. The BH naturally emerges from
the theory of Gravity called as General Relativity [2, 3]. The Kerr metric [4]
describes the geometry of empty space-time around a rotating uncharged axially
symmetric black hole (BH) with a quasi-spherical event horizon. The Kerr
metric is an exact solution of the Einstein field equations of general
relativity. A distinctive prediction of general relativity for the Kerr metric
is that a rotating body should exhibit frame-dragging known as Lense–Thirring
precession. In the case of Kerr BH at close enough distances, all objects even
light must rotate with the BH; the region where this holds is called the
ergosphere. The Kerr BH have surfaces where the metric seems to have apparent
singularities; the size and shape of these surfaces depends on the BH’s mass
and angular momentum. The outer surface encloses the ergosphere and has a
shape similar to a flattened sphere. The inner surface marks the event horizon
– objects passing into the interior of this horizon can never again
communicate with the world outside that horizon. Objects between these two
surfaces must corotate with the rotating BH. This feature can in principle be
used to extract energy, from a rotating BH, equal to its invariant mass
energy.

In some extension of Einstein gravity the Kerr black hole has been studied
recently in [5, 6]. Nearly static magnetized Kerr black-hole in nonlinear
electrodynamics has been studied recently in [7, 8].

In this paper we found the radial equation for equatorial motion for Kerr BH
in terms of the effective potential. We found the solution of the radial
equation of the Kerr metric and found the expression of the angular momentum
per unit mass squared, $\tilde{a}=\frac{a}{M}$. We showed the variations of
$\tilde{a}$ as a function of $\frac{r}{M}$. Finally we enumerated the maximum
values of radius of innermost stable circular orbit ($r_{ISCO}$).

The paper is organized is as follows. In Sec.2 we studied the Kerr geometry
and the effective potential. In Sec.3 we have shown the solutions to the
radial equation of the Kerr Metric. The numerical studies and correlation
among various parameters have been analyzed in Sec.4. We conclude in Sec.5.

##  2 Kerr Geometry

For the Kerr BH with mass $M$, angular momentum $J$ and the
coordinates$(t,r,\theta,\phi)$ the line element of space-time (in geometrized
units with $G=c=1$) is given by[4, 9]:

${\displaystyle
ds^{2}=-(1-\frac{2Mr}{\rho^{2}})dt^{2}-\frac{4Mar\sin^{2}{\theta}}{\phi^{2}}d\phi
dt+\frac{\rho^{2}}{\Delta}dr^{2}+\rho^{2}d\theta^{2}+\newline
(r^{2}+a^{2}+\frac{2Mra^{2}\sin^{2}{\theta}}{\rho^{2}})\sin^{2}{\theta}d\phi^{2}}$
(1)

where,
$\rho^{2}=r^{2}+a^{2}\cos^{2}{\theta}$ ,
$\Delta=r^{2}-2Mr+a^{2}$
$a=\frac{J}{M}$
It is obvious from the line-element that the corresponding tensor is not the
diagonal form due to the presence $d\phi dt$ term. We see that the tensor
corresponding to the line element in Eq.1 is given by:

${\displaystyle
g_{\mu\nu}=\begin{pmatrix}-(1-\frac{2Mr}{\rho^{2}})&0&0&-\frac{2Mar\sin^{2}\theta}{\rho^{2}}\\\
0&\frac{\rho^{2}}{\Delta}&0&0\\\ 0&0&\rho^{2}&0\\\
-\frac{2Mar\sin^{2}\theta}{\rho^{2}}&0&0&(r^{2}+a^{2}+\frac{2Ma^{2}r\sin^{2}\theta}{\rho^{2}})\sin^{2}{\theta}\end{pmatrix}}$
(2)

with off-diagonal elements $g_{t\phi}=g_{\phi t}$. The metric is a solution of
the vacuum Einstein equation and we shall explore some of its most important
properties.

It follows from the Kerr metric that it is independent of $t$ and hence
stationary. The metric is also independent of $\phi$ thus there exist Killing
vectors $\xi^{\mu}=(1,0,0,0)^{\tau}$ and $\eta^{\mu}=(0,0,0,1)^{\tau}$. The
Eq.1 reveals singularities when $\rho$ and $\Delta$ vanishes. The singularity
for $\rho=0$ arises when $r=0$ and $\theta=\frac{\pi}{2}$, which corresponds
to a real, physical singularity and is the generalization of the physical
singularity in the Schwarzschild metric [9].

We have two solutions for $\Delta=0$:

${\displaystyle r\pm=M\pm\left(M^{2}-a^{2}\right)^{\frac{1}{2}}}$ (3)

The event horizon does not exist for $a>M$ (as this is a coordinate
singularity for Schwarzschild metric). Furthermore it is clear that under the
condition $a=0$, Eq.3 reduces to $r_{+}=2M$ (i.e., event horizon in the
Schwarzschild geometry). Thus Eq.3 defines the generalization of the event
horizon. Evidently then the event horizon ranges from $r=2M$ (when $a=0$) to
$r=M$ (when $a=M$ for extremal black holes). It can be saying that as the
angular momentum per mass of black holes increases, the radius of the event
horizon gets smaller. From Eq.3 we also find that the Kerr parameter $a\leq
M$, a limitation which does not apply to other astronomical objects such as
stars.

In general, the orbits of particles in the Kerr geometry are not confined to a
plane unlike in Schwarzschild geometry, due to absence of spherical symmetry.
However the metric is still asymmetric so there exist orbits in the equatorial
plane, and such orbits are especially interesting since the accretion disk
lies close to the equatorial plane[11]. We will restrict our attention to such
orbits and for the equatorial plane we then have $u^{\theta}=0$.

Following Eq.1 the Kerr geometry reduces to:

${\displaystyle ds^{2}=-\left(1-\frac{2M}{r}\right)dt^{2}-\frac{4Ma}{r}d\phi
dt+\frac{r^{2}}{\Delta}dr^{2}+\left(r^{2}+a^{2}+\frac{2Ma^{2}}{r}\right)d\phi^{2}}.$
(4)

Since the Kerr metric is independent of $t$ and $\phi$ there exist two
conserved quantities.

${\displaystyle-e=g_{\mu\nu}\xi^{\mu}u^{\nu}~{}~{}~{}~{}~{}~{}and~{}~{}~{}~{}~{}~{}l=g_{\mu\nu}\eta^{\mu}u^{\nu}}$
(5)

where $e$ is the energy per unit rest mass and $l$ is the angular momentum per
unit rest mass. The $\xi^{\mu}$ and $\eta^{\mu}$ are two Killing vectors. From
the Kerr metric it is clear that both the quantities in Eq.5 are linear
combinations of $u^{t}$ and $u^{\phi}$ according to

${\displaystyle-e=g_{tt}u^{t}+g_{t\phi}u^{\phi}}$ (6) ${\displaystyle
l=g_{\phi t}u^{t}+g_{\phi\phi}u^{\phi}}$ (7)

One can invert to solve $u^{t}$ and $u^{\phi}$ and are as follows:

${\displaystyle u^{t}=-\frac{1}{\Delta}(lg_{t\phi}-eg_{\phi\phi})}$ (8)
${\displaystyle u^{\phi}=-\frac{1}{\Delta}(lg_{tt}+eg_{\phi t})}$ (9)

One can written out more explicitly as:

${\displaystyle\frac{dt}{d\tau}=\frac{1}{\Delta}\left[\left(r^{2}+a^{2}+\frac{2Ma^{2}}{r}\right)e-\frac{2Ma}{r}l\right]}$
(10)
${\displaystyle\frac{d\phi}{d\tau}=\frac{1}{\Delta}\left[\left(1-\frac{2M}{r}\right)l+\frac{2Ma}{r}e\right]}$
(11)

Now applying normalization condition $g_{\mu\nu}u^{\mu}u^{\nu}=-1$ one can
find the radial equation for equatorial motion:

${\displaystyle\frac{e^{2}-1}{2}=\left(\frac{dr}{d\tau}\right)^{2}+V_{eff}(r)}$
(12)

where $V_{eff}(r)$ is the effective potential given by:

${\displaystyle
V_{eff}(r)=\frac{l^{2}-a^{2}(e^{2}-1)}{2r^{2}}-\frac{M(l-ae)^{2}}{r^{3}}-\frac{M}{r}}$
(13)

We’ve plotted the $V_{eff}$ as a function of the radial distances (in
lightyear) in Fig.1 for different values of a. We see that $\frac{r}{M}\geq
15$ the values of the effective potential is insensitive to the rotation
parameter.

Figure 1: In the left-panel $V_{eff}$ as a function of $\frac{r}{M}$ for Kerr
geometry. The blue, red and black are corresponds to the value of
$a=0.1M,0.5M,0.9M$ respectively for $\ell=3M$. In the right panel is the same
for $\ell=2\sqrt{3}M$.

Figure 2: The $V_{eff}$ as a function of $\frac{r}{M}$ and $\frac{\ell}{M}$
(where $r$ is stands for $r_{ISCO}$ (in light years), M is in Solar mass) and
$\ell$ is angular momentum per unit mass.

The three dimensional profile of effective potential is shown in Fig.2 as a
function of $\frac{r}{M}$ and $\frac{\ell}{M}$. The $r$ is stands for
$r_{ISCO}$ (in light years) and M is in Solar mass and $\ell$ is angular
momentum per unit mass.

Figure 3: $V_{eff}$ as a function of $r$ and $a=\frac{\ell}{M}$ for M=100
Solar mass with different orientations.

In Fig.3 we have shown three dimensional effective potential profile as a
function of $r$ and $a=\frac{\ell}{100}$. We showed in different orientations
so that the $v_{eff}$ profile is more clear from different angle.

##  3 Solutions of Radial Equation of the Kerr Metric

In principle we may imagine that values for $r_{ISCO}$ (radius of innermost
stable circular orbit) are known and that the angular momentum per mass
squared $\tilde{a}=\frac{a}{M}$ is an unknown which one wish to solve for.

Disk-accretion affects $\tilde{a}$. As black holes accrete matter and in order
to do so expressions for the angular momentum per unit rest mass of particles
($l$) and the energy per unit rest mass of particles ($e$) are very important.
So, solving Eqn.12 for $e$ and $\tilde{a}$ in terms of $r_{ISCO}$ one can
define the potential $V(r)$ as:

$V(r)=-2r^{2}\left(V_{eff}(r)-\frac{e^{2}-1}{2}\right)$

${\displaystyle V(r)=(r^{2}+a^{2})(e^{2}-1)+\frac{2M(l-ae)^{2}}{r}+2Mr-l^{2}}$
(14)

where $r>0$ and here $r_{ISCO}$ is the value of r for which we get the
innermost stable, or marginally stable orbit with constant $r$. For $r_{ISCO}$
it necessarily follows $\frac{dr}{d\tau}=0$ and for particles to remain in
circular orbit the conditions $V(r)=0$ and $V_{0}(r)=0$ must hold. For a
stable orbit we have the condition $V^{\prime\prime}(r)>0$, and since
$r_{ISCO}$ is the orbit which is on the verge of being unstable one can find
that this condition becomes an equality: $V^{\prime\prime}(r)=0$. The
derivative of $V(r)$ is given by:

${\displaystyle V^{\prime}(r)=2r(e^{2}-1)-\frac{2M(l-ae)^{2}}{r^{2}}+2M}$ (15)

If we take the double derivative of Eqn.15 we get

${\displaystyle V^{\prime\prime}(r)=2(e^{2}-1)+\frac{4M(l-ae)^{2}}{r^{3}}}$
(16)

Now multiply of Eqn.15 with r and add with Eqn.14 we can get much simpler
expression 111Using computer algebraic methods one can solve energy ($e$) and
angular momentum ($l$) in terms of $r$, $M$ and $a$(angular momentum per unit
mass) [10].

${\displaystyle V(r)+rV^{\prime}(r)=(3r^{2}+a^{2})(e^{2}-1)+4Mr-l^{2}}$ (17)

where $r>0$ and here $r_{ISCO}$ is the value of r of the innermost stable
circular orbit(ISCO). As $r_{ISCO}$ is the minimum stable point so
$V^{\prime}(r)=0$ and for the innermost stable circular orbit in the
equatorial plane the second derivative $V^{\prime\prime}(r)=0$, therefore we
can differentiate Eq.17 second time and get

${\displaystyle\frac{d}{dr}[V(r)+rV^{\prime}(r)]=2V^{\prime}(r)+rV^{\prime\prime}(r)}$
(18)

The above equation only contains the first and second derivative.

As soon as we have imposed the condition $V^{\prime\prime}(r)=0$ the radius
$r$ must necessarily be $r_{ISCO}$. Hereafter whenever we use $r$ it always
referring to $r_{ISCO}$.

With this in mind we find that Eq.18 reduces to:

${\displaystyle\frac{d}{dr}[V(r)+rV^{\prime}(r)]=6r(e^{2}-1)+4M=0}$ (19)

and from this it follows that

${\displaystyle e=\left(1-\frac{2M}{3r}\right)^{\frac{1}{2}}}$ (20)

considering only the positive solution. The energy ($e$) per unit rest mass as
a function of $\frac{r}{M}$ is plotted in the left-panel of Fig.4.

In the right-panel of Fig.4 we’ve plotted $\ell$ as a function of $r$ for
M=1,10,100 Solar mass. It is clear from the figure that the finite values of
$\ell$ is started from 0.66, 6.6 and 66.6 for M=1,10,100 respectively shows
that it simply scales.

Figure 4: In the left-panel $e$ (energy per unit rest mass) is plotted as a
function of $\frac{r}{M}$ following Eqn.20. Please note that the values of
$\frac{r}{M}\sim 0.66$ for real $e$. In the right-panel $\ell$ is plotted as a
function of $r$ for M=1(dashed blue), M=10 (dotted-dash red) and M=100(solid
black). The finite value of $\ell$ is starting from 0.66, 6.6 and 66.6 for
M=1,10,100 respectively.

For the critical point Eqn.16 is equal to zero and from Eqn.20 one can get the
expression for $\ell$

${\displaystyle l=a\sqrt{(}1-\frac{2M}{3r})\pm{\frac{r}{\sqrt{3}}}}$ (21)

From Eqn.17 and Eqn.20 and doing simple algebra one can get

${\displaystyle l^{2}=2M(r-\frac{a^{2}}{3r})}$ (22)

Now from Eqn.21 and Eqn.22 with simple algebra one can get second order
equations of $a$:

${\displaystyle
a^{2}\pm{\frac{2r}{\sqrt{3}}}\sqrt{(}1-{\frac{2M}{3r}})a-rM+\frac{r^{2}}{6}=0}$
(23)

The solution of this Eqn.23 is

${\displaystyle
a=\pm\frac{4\sqrt{Mr}}{3}\mp\frac{r}{\sqrt{3}}\sqrt{\left(1-\frac{2M}{3r}\right)}}$
(24)

By defining $\tilde{a}=\frac{a}{M}$ Eq.24 can be simplified as:

${\displaystyle\tilde{a}=\pm\frac{4}{3}\sqrt{\frac{r}{M}}\mp\frac{r}{M\sqrt{3}}\sqrt{1-\frac{2M}{3r}}}$
(25)

where the upper signs refer to corotation and the lower signs refer to counter
rotation.

Figure 5: The three dimensional profile of $a$ as a function of $r=[0:10]$ and
$M=[1-100]$ (in unit of Solar mass) with positive (left-panel) and
negative(right-panel) sign for the first-term of Eqn. .

We’ve shown in Fig.5 the three dimensional profile of $a$ as a function of $r$
and $M$ with positive (left-panel) and negative (right-panel) sign for the
first-term of Eqn.24.

The profile of $\tilde{a}$ as a function of $\frac{r}{M}$ is plotted in the
right-panel of Fig.2.

Figure 6: $\tilde{a}$ as a function of $\frac{r}{M}$ following Eqn.25 in the
left-panel. In the right-panel $\frac{\ell}{M}=\frac{J}{M^{2}}$ is plotted as
a function of $\frac{r}{M}$ following Eqn.26 taking +ve sign only.

We’ve plotted the angular momentum per unit mass square mass, i.e.,
$\tilde{a}=\frac{\ell}{M^{2}}$ as a function of $\frac{r}{M}$ following Eqn.25
in the left-panel of Fig.6. It is clear that a counter rotates particles
starts to corotate with Black hole once $\frac{r}{M}\geq 6$ and vice-versa. We
can say that is the outer most boundary of ergosphere for the upper sign of
$\tilde{a}$.

In the right-panel of Fig.6 we plotted Eqn.26 (i.e., the angular momentum per
unit mass) as a function of $\frac{r}{M}$. Lower than the values of
$\frac{r}{M}\leq 0.66$ is not physically not possible because $M$ is increased
very much so to rotate the black hole needs huge energy say nearly infinity
and is not possible physically.

Figure 7: $V_{eff}$ as a function of $\frac{r}{M}$ and $\tilde{a}=\frac{a}{M}$
for $\frac{\ell}{M}$= 3, 2 $\sqrt{3}$ and 4 respectively for left, middle and
right panel. Please note that the values $V_{eff}$ is zero if $\frac{r}{M}\leq
0.66$ for all values of $\tilde{a}=\frac{a}{M}$.

In Fig.7 we’ve plotted $V_{eff}$ as a function of $\frac{r}{M}$ and
$\tilde{a}=\frac{a}{M}$ for $\frac{\ell}{M}$= 3, 2 $\sqrt{3}$ and 4
respectively for left, middle and right panel. It is clear that the values
$V_{eff}$ is zero if $\frac{r}{M}\leq 0.66$ for all values of
$\tilde{a}=\frac{a}{M}$. We see that for a particular value of $\tilde{a}$ the
value of $V_{eff}$ create the barrier in small region of $\frac{r}{M}$ which
prevent the object to collapse into black hole. The second term of $V_{eff}$
is dominant with respect to the first term for the low values of
$\frac{r}{M}$.

##  4 Numerical Analysis

Finally we expresses $\ell$ as a function of $r$ using the expressions of
$\tilde{a}$ and $e$. Using Eq.21 we solved for $\ell$. Taking $a=\tilde{a}M$
and the upper sign and with simplification leads to

${\displaystyle
l=\frac{2M}{3^{\frac{3}{2}}}\left[1+2\left(\frac{3r}{M}-2\right)^{\frac{1}{2}}\right]}$
(26)

We have now derived three expressions for $e$, $l$ and $\tilde{a}$ in terms of
$z$, where $z=\frac{r_{ISCO}}{M}$ is called the Boyer-Lindquist radius [12].

${\displaystyle e=\left(1-\frac{2}{3z}\right)^{\frac{1}{2}}}$ (27)

${\displaystyle
l(z.M)=\frac{2M}{3^{\frac{3}{2}}}\left[1+2\left(3z-2\right)^{\frac{1}{2}}\right]}$
(28)

${\displaystyle\tilde{a}=\pm\frac{1}{3}\sqrt{z}\left[4-\sqrt{\left(3z-2\right)}\right]}$
(29)

We mention here that $z$ can be solved for as a function of $\tilde{a}$
following [10]:

${\displaystyle
z(\tilde{a})=3+w_{2}\mp\left[(3-w_{1})(3+w_{1}+2w_{2})\right]^{\frac{1}{2}}}$
(30)

where

${\displaystyle
w_{1}=1+(1-\tilde{a}^{2})^{\frac{1}{3}}\left[(1+\tilde{a})^{\frac{1}{3}}+(1-\tilde{a})^{\frac{1}{3}}\right]}$
(31)

${\displaystyle w_{2}=(3\tilde{a}+w_{1}^{2})^{\frac{1}{3}}}$ (32)

The upper sign in Eqn.30 refers to corotation and the lower sign refers to
counter rotation.

Figure 8: $\tilde{a}$ as a function of $w_{1}$ following Eqn. in the left-
panel. In the right-panel $\tilde{a}$ as a function of $w_{2}$ following Eqn.

In Fig.8 we’ve shown the variations of $\tilde{a}$ as a function of $w_{1}$
and $w_{2}$ in the left-panel and right-panel respectively. This shows the
maximal allowed values of $w_{1}$ and $w_{2}$ in the Schwarzschild limit. This
limiting values of $w_{1}$ and $w_{2}$ have been taken into considerations
while plotting the three dimensional profile of $z(\tilde{a})$, i.e., the
Boyer-Lindquist radius in Fig.10.

Figure 9: $z(\tilde{a})$ as a function of $\tilde{a}$ following Eqn.29 with
negative (left-panel) and positive(right-panel) signature.

We’ve shown the Boyer-Lindquist radius ($z$), i.e., the radius of rotating
particle around the black hole, as a function of $\tilde{a}$ following Eqn.29
with negative and positive sign convention in the left and right-panel of
Fig.9 respectively considering $M=1$. The Boyer-Lindquist radius is maximum
for $\tilde{a}=0$ and this is equivalent to the Schwarzschild case. The space
time fabric is directly converges towards black hole – so energy required to
rotate around black hole is the kinetic energy of the particle itself. Now
once $\tilde{a}$ increase the space time fabric get twisted so as a result the
geodesic path of particle also twisted around the BH. In this scenario the
particle need less amount of kinetic energy for the stable orbit as compared
to Schwarzschild case. For $\tilde{a}$ =1 the maximum value of innermost
stable orbit nearly $z\simeq 8.4(1.0)$ for positive and negative signature in
Eqn.29.

Figure 10: The three dimensional profile of $z(\tilde{a})$ as a function of
$w_{1}$ and $w_{2}$ with positive (left-panel) and negative(right-panel) sign
convention following Eqn… The maximal values of $w_{1}$ and $w_{2}$ are follow
from Fig.8.

In Fig.10 we shown the three dimensional profile of $z(\tilde{a})$ with
$w_{1}$ and $w_{2}$ with positive (left-panel) and negative (right-panel) sign
convention following Eqn.30. The allowed maximal values of $w_{1}$ and $w_{2}$
are taken from Fig.8.

##  5 Summary and Conclusion

The BH naturally emerges from the metric solution of the general relativity.
The rotating BH as known as Kerr and general relativity put an upper limit on
the angular momentum per mass squared of black holes $\leq 1$. If the value is
more than one then above which the event horizon of the Kerr BH is not exist.
In this paper we have studied the effective potential of the rotating Kerr
metric and shown the effective potential profile for different values of
rotation parameter ($a$). The solution of the radial equation of the Kerr
metric is found and hence of the angular momentum per unit mass squared
$\tilde{a}$. The behaviors of $\tilde{a}$, the energy per unit rest mass ($e$)
are shown as a function of $\frac{r}{M}$. We demarcated the maximum values of
radius of innermost stable circular orbit, i.e., $r_{ISCO}$, for
$\tilde{a}=1$. Our numerical findings are relevant while studying the BH’s
accretion of the matter from surrounding areas, consequently increasing the
mass and the angular momentum of the BHs. The solutions of the angular
momentum equations are very important to study the accreation of the matter
into the BHs and we will address this issues in other study.

## References

  * [1] S. Vagnozzi, R. Roy, Y. D. Tsai, L. Visinelli, M. Afrin, A. Allahyari, P. Bambhaniya, D. Dey, S. G. Ghosh and P. S. Joshi, et al. Class. Quant. Grav. 40, no.16, 165007 (2023) doi:10.1088/1361-6382/acd97b [arXiv:2205.07787 [gr-qc]].
  * [2] James B. Hartle, Gravity: An Introduction to Einstein’s General Relativity, Addison Wesley, ISBN: 0–8053–8662–9, 2003.
  * [3] Steven Weinberg, Gravitation and Cosmology: Principles and Applications of the Gen- eral Theory of Relativity, John Wiley & Sons, ISBN: 0–471–92567–5, 1972.
  * [4] Roy Patrick Kerr, Gravitational Field of a Spinning Mass as an Example of Alge- braically Special Metrics, Phys. Rev. Letters, Vol. 11, p. 237, 1963.
  * [5] O. Donmez, [arXiv:2307.11725 [astro-ph.HE]].
  * [6] O. Donmez, F. Dogan and T. Sahin, Universe 8, no.9, 458 (2022) doi:10.3390/universe8090458 [arXiv:2205.14382 [astro-ph.HE]].
  * [7] K. G. Managave, H. A. Redekar, R. B. Kumbhar, S. P. Das and K. Y. Rajpure, Int. J. Mod. Phys. A 38, no.28, 2350152 (2023) doi:10.1142/S0217751X2350152X [arXiv:2303.07736 [gr-qc]].
  * [8] H. A. Redekar, R. B. Kumbhar, S. P. Das and K. Y. Rajpure, [arXiv:2308.12639 [gr-qc]].
  * [9] James M. Bardeen, Kerr Metric Black Holes, Nature, Vol. 226, p. 64, 1970.
  * [10] James M. Bardeen, William H. Press and Saul A. Teukolsky, Rotating Black Holes: Locally Nonrotating Frames, Energy Extraction, and Scalar Synchrotron Radiation, The Astrophysical Journal, Vol. 178, p. 347, 1972.
  * [11] Brandon Carter, Global Structure of the Kerr Family of Gravitational Fields, Physical Review, Vol. 174 5, p. 1559–1571, 1968
  * [12] Robert H. Boyer & Richard W. Lindquist, Maximal Analytic Extension of the Kerr metric, Journal of Mathematical Physics, Vol. 8, p. 265, 1967.