# Efficient Reinforcement Learning for Jumping Monopods

Riccardo Bussola1, Michele Focchi1, Andrea Del Prete2, Daniele Fontanelli2,
Luigi Palopoli1 1 The authors are with the Dipartimento di Ingegneria and
Scienza dell’Informazione (DISI), University of Trento<EMAIL_ADDRESS>The authors are with the Dipartimento di Ingegneria Industriale (DII),
University of Trento<EMAIL_ADDRESS>
###### Abstract

In this work, we consider the complex control problem of making a monopod
reach a target with a jump. The monopod can jump in any direction and the
terrain underneath its foot can be uneven. This is a template of a much larger
class of problems, which are extremely challenging and computationally
expensive to solve using standard optimisation-based techniques. Reinforcement
Learning (RL) could be an interesting alternative, but the application of an
end-to-end approach in which the controller must learn everything from
scratch, is impractical. The solution advocated in this paper is to guide the
learning process within an RL framework by injecting physical knowledge. This
expedient brings to widespread benefits, such as a drastic reduction of the
learning time, and the ability to learn and compensate for possible errors in
the low-level controller executing the motion. We demonstrate the advantage of
our approach with respect to both optimization-based and end-to-end RL
approaches.

##  I Introduction

Legged robots have become a popular technology to navigate unstructured
terrains, but are complex devices for which control design is not trivial.
Remarkable results have been reached for locomotion tasks like walking and
trotting [1]. Other tasks, like performing jumps, are more challenging because
even small deviations from the desired trajectory can have a large impact on
the landing location and orientation [2]. This problem has received some
attention in the last few years. A line of research has produced heuristic
approaches relying on physical intuitions and/or on simplified models to be
used in the design of controllers or planners [3, 4]. However, the hand-
crafted motion actions produced by these approaches are not guaranteed to be
physically implementable. Another common approach is to use full-body
numerical optimisation [5]. Very remarkable is the rich set of aerial motions
produced by MIT Mini Cheetah in [6, 7, 8] using a centroidal momentum-based
nonlinear optimisation. A problem with optimisation-based approaches on high
dimensional nonlinear problems is the high computational cost, which makes
them unsuitable for a real-time implementation, especially to replan
trajectories over a receding horizon. Recent advances [9, 5] have brought to
significant improvements in the efficiency of Model Predictive Control (MPC)
for jumping tasks. However, the price to pay is to introduce some artificial
constraints such as fixing the contact sequence, the time-of-flight, or
optimising the contact timings offline.

A third set of approaches is based on Reinforcement Learning (RL). The seminal
work of Lillicrap [10] showed that a Deep Deterministic Policy Gradient (DDPG)
algorithm combined with a Deep Q network could be successfully applied to
learn end-to-end policies for a continuous action domain, using an actor-
critic setting. In view of these results, several groups have then applied RL
to quadrupeds for locomotion tasks [11, 12, 13, 14, 15], and to in-place
hopping legs [16]. As with most of model-free reinforcement approaches, DDPG
requires a large number of training steps (on the order of millions) to find
good solutions. Other approaches [17, 18, 19] seek to improve the efficiency
and the robustness of the learning process by combining Trajectory
Optimization (TO) with RL: they use the former to generate initial
trajectories to bootstrap the exploration of alternatives made by the latter.

As a final remark, the efficiency and the robustness of the RL learning
process is heavily affected by a correct choice of the action space [20, 21].
Some approaches require that the controller directly generates the torques
[22], while others suggest that the controller should operate in Cartesian or
joint space [12, 23].

Paper Contribution. This work proposes an RL framework capable of producing an
omni-directional jump trajectory (from standstill) on uneven terrain, computed
within a few milliseconds, unlocking real-time usage with the current
controller rates (e.g. in the order of $kHz$). The main objective of this
paper is to reduce the duration of the learning phase without sacrificing the
system’s performance. A reduced length of the learning phase has at least two
indisputable advantages: lowering the barriers for the access to this
technology for professionals and companies with a limited availability of
computing power, and addressing the environmental concerns connected with the
carbon footprint of learning technologies [24].

Our strategy is based on the following ideas: first, learning is performed in
Cartesian space rather than in joint space so that the agent can more directly
verify the effect of its actions. Second, we notice that while the system is
airborne, its final landing point are dictated by simple mechanical laws
(ballistic). Therefore, the learning process can focus solely on the thrusting
phase. Third, we know from biology [25] that mammals are extremely effective
in learning how to walk because of ”prior” knowledge in their genetic
background. This means that the learning process can be guided by an
approximate knowledge of what the resulting motion should ”look like”.
Specifically, we parametrise the thrusting trajectory (i.e. from standstill to
lift-off) for the Center of Mass (CoM) with a ($3^{rd}$ order) Bezier curve.
This choice of a Bezier curve to parametrize the action space is not uncommon
in the literature [26], [27], [28], and in our case is motivated by the
physical intuition that a jump motion, to exploit the full joint range, needs
a ”charging” phase to compress the legs followed by an extension phase where
the CoM is accelerated both upwards and in the jump direction. Clearly, by
making this restriction, we prevent the learning phase from exploring
alternative options. However, as we will see in the paper, the final result is
very close to optimal and the system retains good generalisation abilities.

Our approach is based on TD3, a state-of-the-art Deep Reinforcement Learning
(Deep-RL) algorithm [29] trained to minimise a cost very similar to the ones
typically used in optimal control. The Cartesian trajectory generated by our
RL agent is translated into joint space via inverse kinematics, and tracked by
a low-level joint-space Proportional-Derivative (PD) controller with gravity
compensation. Our results reveal that possible inaccuracies in the controller
can be learned and compensated by the RL agent. In this paper we do not focus
on the landing phase, which we assume to be managed by a different controller,
such as [2]. We compare our approach (that we will call Guided Reinforcement
Learning (GRL)) with both a baseline TO controller with a fixed duration of
the thrusting phase, and a ”standard” End-to-end Reinforcement Learning (E2E)
(which considers joint references as action space). In the first case, we
achieved better or equal performance, and a dramatic reduction in the online
computation time. With respect to E2E locomotion approaches [22, 12], we
observed a substantial reduction in the number of episodes (without
considering parallelization) needed to achieve good learning performances.
Instead, an E2E implementation specific for a jump motion (Section IV-2), did
not provide satisfactory learning results.

The paper is organized as follows: Section II presents the GRL approach,
detailing the core components of the MDP (action, reward functions). Section
III provides implementation details. In Section IV we showcase our simulation
results comparing with state-of-the-art approaches. Finally, we draw the
conclusions in Section V.

##  II Problem description and solution overview

Simple notions of bio-mechanics suggest that legged animals execute their
jumps in three phases: 1\. _thrust:_ an initial compression is followed by an
explosive extension of the limbs in order to gain sufficient momentum for the
lift-off; the phase finishes when the foot leaves the ground; 2\. _flight:_
the body, subject uniquely to gravity, reaches an apex where the vertical CoM
velocity changes its sign and the robot adjusts its posture to prepare for
landing; 3\. _landing:_ the body realizes a touch-down, which means that the
foot establishes again contact with the ground. For the sake of simplicity, we
consider a simplified and yet realistic setting: a monopod robot, whose base
link is sustained by passive prismatic joints preventing any change in its
orientation (see Section III-C). The extension to a full quadruped (which
requires considering also angular motions) is left for future work. In this
paper, _we focus on the thrust phase_. The flight phase is governed by the
ballistic law. Let $\mathbf{c}_{tg}$ be the target location and let the CoM
state at lift-off be $(\mathbf{c}_{lo},\mathbf{\dot{c}}_{lo})$. After lift-
off, the trajectory lies on the vertical plane containing $\mathbf{c}_{lo}$
and $\mathbf{c}_{tg}$. The set of possible landing CoM positions are function
of $\mathbf{c}_{lo}$, $\mathbf{\dot{c}}_{lo}$ and are given by the following
equation:

$\begin{cases}\mathbf{c}_{tg,xy}=\mathbf{c}_{lo,xy}+\mathbf{\dot{c}}_{lo,xy}T_{fl}\\\
\mathbf{c}_{tg,z}=\mathbf{c}_{lo,z}+\mathbf{\dot{c}}_{lo,z}T_{fl}+-\frac{1}{2}gT_{fl}^{2}\end{cases}$
(1)

where
$T_{fl}=(\mathbf{c}_{tg,xy}-\mathbf{c}_{lo,xy})/\mathbf{\dot{c}}_{lo,xy}$ is
the flight time. In this setting, our problem can be stated as follows:

###### Problem 1

Synthesise a thrust phase that produces a lift-off configuration (i.e. CoM
position and velocity) that: 1. satisfies (1), 2\. copes with the potentially
adverse conditions posed by the environment (i.e. contact stability, friction
constraints), 3\. satisfies the physical and control constraints.

Nonlinear optimisation is frequently used for similar problems. However, it
has two important limitations that discourage its application in our specific
case: 1. the computation requirements are very high, complicating both the
real–time execution and the use of low–cost embedded hardware, 2. the problem
is strongly non-convex, which can lead the solvers to local minima.

###  II-A Overview of the approach

In this work, we use RL to learn optimal joint trajectories to realise a jump
motion, that is then tracked by a lower-level controller. We adopted the
state-of-the-art Twin Delayed Deep Deterministic Policy Gradient (TD3), a
Deep-RL technique. This algorithm is based on the well-known Actor-Critic
architecture. The main idea behind TD3 is to adopt two Deep Neural Network
(NN) to approximate the policy, represented by the Actor, and two Deep NN to
approximate the action-value function, represented by the Critic.

Figure 1: Diagram of the RL Framework. The framework is split into two levels:
the RL agents and the planner. The RL agent produces an action for the
planner, based on a desired target. This computes a Bezier reference curve
that is mapped into joint motion via inverse kinematics and tracked by the PD
controller that provides the joint torques to feed the robot. During the
training, at the end of each episode a reward is computed and fed back to the
RL agent.

Our RL pipeline is depicted in Fig. 1.

##  III Design and implementation of the RL agent

Based on the discussion in the previous section, the thrust phase is
characterised by the lift-off position $\mathbf{c}_{lo}$ and velocity
$\mathbf{\dot{c}}_{lo}$ and by the thrust time $T_{th}$ which is the time
spent to reach the lift-off configuration from the initial state. The state of
the environment is defined as ($\mathbf{c}$, $\mathbf{c}_{tg}$) where
$\mathbf{c}\in\mathbb{R}^{3}$ is the CoM position and
$\mathbf{c}_{tg}\in\mathbb{R}^{3}$ the CoM at the landing location (target).
The objective of the RL agent is to find the jump parameters
($\mathbf{c}_{lo}$, $\mathbf{\dot{c}}_{lo}$, $T_{th}\in\mathbb{R}$), that
minimise the landing error at touch-down $\|\mathbf{c}-\mathbf{c}_{tg}\|$
while satisfying the physical constraints. Our jumping scenario thought can be
seen as a single-step trajectory where the only action performed leads always
to the end state.

###  III-A The Action Space

The dimension of the action space has a strong impact on the performance of
the RL algorithm. Indeed, a NN with a smaller number of outputs is usually
faster to train. What is more, a smaller action space reduces the complexity
of the mapping, speeding up the learning process.

A first way to reduce the complexity of the action space is by expressing
$\mathbf{c}_{lo}$ and $\mathbf{\dot{c}}_{lo}$ in spherical coordinates.
Because of the peculiar nature of a jump task, the trajectory lies in the
plane containing the CoM $\mathbf{c}$ and its desired target location
$\mathbf{c}_{tg}$. Hence, the yaw angle $\varphi$ remains constant
($\bar{\varphi}$) throughout the flight and we can further restrict the
coordinates to a convex bi-dimensional space:

$\displaystyle\small\begin{cases}\mathbf{c}_{lo,x}=r~{}\cos(\theta)\cos(\bar{\varphi})\\\
\mathbf{c}_{lo,y}=r~{}\cos(\theta)\sin(\bar{\varphi})\\\
\mathbf{c}_{lo,z}=r~{}\sin(\theta)\end{cases}\begin{cases}\mathbf{\dot{c}}_{lo,x}=r_{v}~{}\cos(\theta_{v})\cos(\bar{\varphi})\\\
\mathbf{\dot{c}}_{lo,y}=r_{v}~{}\cos(\theta_{v})\sin(\bar{\varphi})\\\
\mathbf{\dot{c}}_{lo,z}=r~{}\sin(\theta_{v})\end{cases}$ (2)

As shown in Fig. 2, the lift-off position $\mathbf{c}_{lo}$ is identified by:
the radius $r$ (i.e., the maximum leg extension), the yaw angle $\varphi$ and
the pitch angle $\theta$. Likewise, the lift-off velocity $\dot{c}_{lo}$, is
described by its magnitude $r_{v}$, and the pitch angle $\theta_{v}$ with
respect to the ground.

Figure 2: Action parametrization and its bounds. On the left the top view, on
the right the side view of the jumping plane.

Therefore, by using this assumption, we have reduced the dimension of the
action space from $7$ to $5$:
$\mathbf{a}=(T_{th},r,\theta,r_{v},\theta_{v})\in\mathbb{R}^{5}$.

The action space can be further restricted by applying some domain knowledge.
The radius $r$ has to be smaller than a value $r_{max}$ ($0.32$ $m$) to
prevent boundary singularity due to over-extension, and greater than a value
$r_{min}$($0.25$ $m$) to avoid complete leg retraction. The bounds on the
velocity $\dot{c}_{lo}$, represented by $r_{v}\in\left[0.1,4\right]$ $m/s$,
and $\theta_{v}\in\left[\frac{\pi}{6},\frac{\pi}{2}\right]$ $rad$, and the
bounds on pitch angle $\theta\in\left[\frac{\pi}{4},\frac{\pi}{2}\right]$
$rad$ are set to rule out jumps that involve excessive foot slippage and
useless force effort. Specifically, restricting to a positive $\theta_{v,min}$
ensures a non-negligible vertical component for the velocity, while bounding
$\theta_{v}$ to the positive quadrant secures that the lift-off velocity be
oriented ”toward” the target.

####  III-A1 Trajectory Parametrisation in Cartesian Space

Our strategy to tackle the problem of generating a compression-extension
trajectory for the leg, to achieve a given lift-off configuration
$\mathbf{c}_{lo}$, is based on two important choices: 1\. making the RL agent
learn the trajectory in Cartesian Space and then finding the joint
trajectories through inverse kinematics, 2. restricting the search of the
Cartesian space evolution to curves generated by known parametric functions.
This has the aim to reduce the search space and simplifying convergence.

The analytical and geometric properties of $3^{rd}$ order Beziér curves make
them a perfect fit for our problem. A $3^{rd}$ order Bèzier curve is defined
by four control points. In our case, the first and the final points are
constrained to be the initial and final CoM positions, respectively. The
derivative of a $3^{rd}$ degree Bèzier curve is itself a Bèzier curve of
$2^{nd}$ degree with 3 control points defined as
$3(\mathbf{P}_{i+1}-\mathbf{P}_{i})$. The curve domain is defined only in the
normalised time interval: $t\in\begin{bmatrix}0,1\end{bmatrix}$. Defining the
following Bernstein polynomials:

$\displaystyle\bm{\eta}(t)=\begin{bmatrix}(1-t)^{3}&3(1-t)^{2}t&3(1-t)t^{2}&t^{3}\end{bmatrix}^{T}$
(3)
$\displaystyle\bm{\dot{\eta}}(t)=\begin{bmatrix}(1-t)^{2}&2(1-t)t&t^{2}\end{bmatrix}^{T}$
(4)

we can compactly write the Beziér curve as function of its
$\mathbf{P}_{i}\in\mathbb{R}^{3}$ control points:

$\displaystyle\mathbf{c}=\begin{bmatrix}\mathbf{P}_{0}&\mathbf{P}_{1}&\mathbf{P}_{2}&\mathbf{P}_{3}\end{bmatrix}\bm{\eta}(t)$
(5)

Since we are considering an execution time
$T_{exe}\in\begin{bmatrix}0,T_{th}\end{bmatrix}$ and
$t=\frac{T_{exe}}{T_{th}}$, then the derivative writes:

$\displaystyle\mathbf{\dot{c}}=\frac{1}{T_{th}}\begin{bmatrix}\mathbf{P}^{\prime}_{0}&\mathbf{P}^{\prime}_{1}&\mathbf{P}^{\prime}_{2}\end{bmatrix}\bm{\dot{\eta}}(t)$
(6)

From the definition of the curve (5) and its derivative (6), we can compute
the control points $\mathbf{P}_{i}$ by setting the boundary conditions of the
initial/lift-off CoM position $\mathbf{c}_{0}$, $\mathbf{c}_{lo}$ and
initial/lift-off CoM velocity $\mathbf{\dot{c}}_{0}$,$\mathbf{\dot{c}}_{lo}$
in (7).

$\small\begin{cases}\begin{aligned}
{\mathbf{P}}^{\prime}_{0}&=\frac{3}{T_{th}}(\mathbf{P}_{1}-\mathbf{P}_{0})\\\
&=\mathbf{\dot{c}}_{0}=\mathbf{0}\\\
{\mathbf{P}}^{\prime}_{1}&=\frac{3}{T_{th}}(\mathbf{P}_{2}-\mathbf{P}_{1})\\\
{\mathbf{P}}^{\prime}_{2}&=\frac{3}{T_{th}}(\mathbf{P}_{3}-\mathbf{P}_{2})\\\
&=\mathbf{\dot{c}}_{lo}\end{aligned}\end{cases}\begin{cases}\begin{aligned}
\mathbf{P}_{0}&=\mathbf{c}_{0}\\\
\mathbf{P}_{1}&=\frac{T_{th}}{3}{\mathbf{P}}^{\prime}_{0}+\mathbf{P}_{0}\\\
&=\frac{T_{th}}{3}\mathbf{\dot{c}}_{0}+\mathbf{c}_{0}\\\
\mathbf{P}_{2}&=-\frac{T_{th}}{3}{\mathbf{P}}^{\prime}_{2}+\mathbf{P}_{3}\\\
&=-\frac{T_{th}}{3}\mathbf{\dot{c}}_{lo}+\mathbf{c}_{lo}\\\
\mathbf{P}_{3}&=\mathbf{c}_{lo}\\\ \end{aligned}\end{cases}$ (7)

###  III-B A physically informative reward function

In RL, an appropriate choice of the reward function is key to the final
outcome. Furthermore, we can use the reward function as a means to inject
prior knowledge into the learning process. In our case, the reward function
was designed to penalise the violations of the physical constraints while
giving a positive reward to the executions that make the robot land in
proximity of the target point. The constraints that must be enforced
throughout the whole thrust phase are called path constraints. To transform
the violations into costs, we introduce a linear activation function
$A(x,\underline{x},\bar{x})$ of the evaluated constraint, as function of the
value $x$ and its lower and upper limits $\underline{x}$, $\bar{x}$:

$A(x,\underline{x},\bar{x})=\left|\min(x-\underline{x},0)+\max(x-\bar{x},0)\right|$

The output of the activation function is zero if the value is in the allowed
range, the violation otherwise.

Physical feasibility check: Before starting each episode, we perform a sanity
check on the action $\mathbf{a}$: if the vertical velocity is not sufficient
to reach the target height, we abort the simulation returning a high penalty
cost $C_{ph}$. This can be computed by obtaining the time to reach the apex
$T_{fup}=\mathbf{\dot{c}}_{{lo},z}/g$ and substituting it in the ballistic
equation:

$\displaystyle\bar{\mathbf{c}}_{z}(T_{fup})$
$\displaystyle=\mathbf{c}_{{lo},z}+\mathbf{\dot{c}}_{{lo},z}T_{fup}+\frac{1}{2}(-g)T_{fup}^{2}$
(8)

this results in
$\bar{\mathbf{c}}_{z}(T_{fup})=\mathbf{c}_{{lo},z}+\frac{1}{2}\frac{\mathbf{\dot{c}}_{{lo},z}^{2}}{g}$
which is the apex elevation. If
$\mathbf{c}_{tg,z}>\bar{\mathbf{c}}_{z}(T_{fup})$, the episode is aborted.

Unilaterality constraint: in a legged robot, a leg can only push on the ground
and not pull. This is because the component of the force $\mathbf{F}$ along
the contact normal ($Z$ for flat terrains), must be positive.
Friction constraint: To avoid slippage, the tangential component of the
contact force $\left\|\mathbf{F}_{x,y}\right\|$ is bounded by the foot-terrain
friction coefficient $\mu$:
$\left\|\mathbf{F}_{x,y}\right\|\leq\mu\mathbf{F}_{z}$.
Joint range and torque constraints: the three joint $\bar{q}_{i}$ kinematic
limits must not be exceeded. Similarly, each of the joint actuator torque
$\bar{\tau}_{i}$ limits must be respected.
Singularity constraint: the singularity constraint avoids the leg being
completely stretched. During the thrusting phase the CoM $\mathbf{c}$ must
stay in the hemisphere of radius equal to the maximum leg extension. This
condition prevents the robot from getting close to points with reduced
mobility that produce high joint velocities in the inverse kinematic
computation. Even though this constraint is enforced by construction in the
action generation, the actual trajectory might still violate it due to
tracking inaccuracies. If a singular configuration is reached, the episode is
interrupted and a high cost is returned. The costs caused by the violation of
path constraints are evaluated for each time step of the thrust phase and
accumulated into the feasibility cost $C_{f}$. In addition to these path
constraints, we also want to take into account the error between actual and
desired lift-off state. This penalty $C_{lo}$ encourages lift-off
configurations that are better tracked by the controller. Another penalty
$C_{td}$ is introduced when an episode does not result in a touchdown. This is
done to force in-place jumps, and avoid the robot to stay stationary. The
positive component of the reward function is the output of a nonlinear landing
target reward function, which evaluates how close the CoM arrived to the
desired target. This reward grows exponentially when this distance approaches
zero:

$R_{lt}(\mathbf{c},\mathbf{c}_{tg})=\frac{\beta}{k\left\|\mathbf{c}-\mathbf{c}_{tg}\right\|+\epsilon},$
(9)

where $k$ is a gain to encourage jumps closer to the target position, and
$\beta$ is an adjustable parameter to bound the max value of $R_{lt}$ and
scale it. An infinitesimal value $\epsilon$ is added at the denominator to
avoid division by zero. Hence, the total reward function is:

$R=1_{\mathbb{R}^{+}}\left[R_{lt}(\mathbf{c},\mathbf{c}_{tg})-\sum_{i=0}^{n_{c}}C_{i}\right]$
(10)

with $n_{c}=8$, and where $C_{i}$ are the previously introduced feasibility
costs. We decided to perform reward shaping[30], by clamping the total reward
to $\mathbb{R}^{+}$ by mean of an indicator function $1_{\mathbb{R}^{+}}$.
This aims to promote the actions that induce constraint satisfaction.

###  III-C Implementation details

The training of the RL agent and the sim-to-sim validation of the learned
policy was performed on top of a Gazebo simulator. Because we are considering
only translational motions, we modelled a 3 Degrees of Freedom (DoFs) monopod
with three passive prismatic joints attached at the base. These prismatic
joints constrain the robot base’s movements to planes parallel to the ground.
For the sake of simplicity, we also considered the landing phase under the
responsibility of a different controller (e.g., see [2]). Our interest was
simply on the touch-down event, which is checked by verifying that the contact
force exceeds a positive threshold $f_{th}$. Therefore, the termination of the
episode is determined by the occurrence of three possible conditions:
execution timeout, singularity, or touch-down event.

The control policy (default NN) is implemented as a neural network that takes
the state as an input, and outputs the actions. The NN is a multi-layer
perceptron with three hidden layers of sizes 256, 512 and 256 with ReLU
activations between each layer, and with tanh output activation to map the
output between -1 and 1. A low-level PD plus gravity compensation controller
generates the torques that are sent to the Gazebo simulator at 1 kHz. The
joint reference positions at the lift-off are reset to the initial
configuration $\mathbf{q}_{0}$ to enable the natural retraction of the leg and
avoid stumbling. Landing locations at different heights are achieved by making
a 5x5 cm platform appear at the desired landing location only at the apex
moment. 111Making the platform appear only at apex is needed for purely
vertical jumps, because it avoids impacts with the platform during the
trusting phase. The impact of the foot with the platform determines the touch-
down moment and the consequent termination of the episode. To train the RL
agent, the interaction with the simulation environment is needed. The
communication between the planner component and the Gazebo simulator is
managed by the Locosim framework [31]. To interact with the planner, and
consequentially, with the environment, we developed a ROS node called
JumplegAgent, where we implemented the RL agent. The code is available at
222Source code available at https://github.com/mfocchi/jump_rl. During the
initial stage of the training process, the action is randomly generated to
allow for an initial broad exploration of the action space for $N_{exp}$
episodes.

##  IV Simulation Results

In this section, we discuss some simulation results that show the validity of
the proposed approach and compare it with state-of-the-art approaches.

Figure 3: Top-view of the feasibility region: (1-5) for different number of
episodes of the training phase (the number of reachable points is computed for
each $X$,$Y$ pair) and (right) in the case of the baseline FDDP.

We used a computer with the following hardware specifications: CPU AMD Ryzen 5
3600, GPU Nvidia GTX1650 4GB, RAM 16 GB DDR4 3200 MHz. During training we
generated targets inside a predefined training region. This samples are
generated randomly inside a cylinder centered on the robot’s initial CoM, with
a radius from 0 to $0.65$ m and a height from $0.25$ m to $0.5$ m. The size
was selected to push the system to its performance limits. The parameters of
the robot, controller and simulation are presented in Table I.

TABLE I: Controller, Planner, Simulator parameters Variable |  Name | Range
---|---|---
$m$ |  Robot Mass [$kg$] | 1.5
P |  Proportional gain | 10
D |  Derivative gain | 0.2
$\mathbf{q}_{0}$ |  Nominal configuration |  $\begin{bmatrix}0&-0.75&1.5\end{bmatrix}$ [$rad$]
$dT$ |  Simulator time step [$s$] | 0.001
$\tau_{max}$ |  Max torque [$Nm$] | 8
$f_{th}$ |  Touch-down force th. [$N$] | 1
$N_{exp}$ |  Num. of expl. steps | 1280 (GRL), 10e4 (E2E)
$b_{s}$ |  Batch size | 256 (GRL), 512 (E2E)
$n_{exp}$ |  Expl. noise | 0.4 (GRL), 0.3 (E2E)
$tg_{rep}$ |  Landing target repetition | 5 (GRL), 20 (E2E)
$N_{train}$ |  Training step interval | 1 (GRL), 100 (E2E)

####  IV-1 Nonlinear Trajectory Optimisation

The first approach to compare with is a standard optimal control strategy
based on FDDP. FDDP is one of the most efficient optimisation algorithms for
whole-body control [32] because it takes full advantage of the intrinsic
sparse structure of the optimal control problem. The FDDP solver is
implemented with the optimal control library Crocoddyl [33] and uses the
library Pinocchio [34] to enable fast computations of costs, dynamics, and
their derivatives.

For the problem at hand, we discretised the trajectory into $N$ successive
knots with a timestep $dT=0.001$ $s$ to ensure high precision. As decision
variables we chose the joint torques. We split the jump in three phases:
thrusting, flying, and landed. The constraints in FDDP are encoded as soft
penalties in the cost. We encoded friction cones, and tracking of foot
locations and velocities at the thrusting/landed stages, respectively. We
added a tracking cost for the CoM reference during the flying phase to
encourage the robot to lift-off. We regularized control inputs and states
throughout the horizon.

####  IV-2 End-to-end RL

At the opposite side of the spectrum of optimal control is the application of
approaches entirely based on Deep Learning, i.e., using RL end-to-end without
injecting any prior domain knowledge. The RL agent sets joint position
references to a low-level (PD)+gravity compensation controller. The use of a
PD controller allows the system to inherit the stabilisation properties of the
feedback controller, but at the same time, it allows for explosive torques (by
regulating the references to have a bigger error w.r.t. to the actual
positions). We query the action until the apex moment because we need to have
set-points for the joints also when the leg is air-borne. After apex, we set
the default configuration $\mathbf{q}_{0}$ for the landing. To be more
specific, instead of directly setting the joint references $\mathbf{q}^{d}$,
the control policy produces as action $\mathbf{a}$ joint angle deviations
$\mathbf{\tilde{q}}\in\mathbb{R}^{3}$ w.r.t. to the nominal joint angle
configuration $\mathbf{q}_{0}$.

As suggested by [35, 36] to ease the learning, we run at a frequency that is
1/5th than the controller one (1 $kHz$). We aggregate the reward for each
control loop iteration and perform a training every $N_{train}=100$ queries.
We terminate the episode if: 1) a touchdown is detected, 2) the robot has
fallen (i.e. the base link close to the ground), 3) we reach singularity 4) a
timeout of 2.5 $s$ is reached. We include in the state CoM and joint
position/velocity. Since the domain is not changing (and the state is
Markovian), augmenting the state with the history of some past samples [37,
12] was not necessary, therefore we tried to keep the problem dimensionality
as low as possible. Hence, the observation state becomes ($\mathbf{c}$,
$\mathbf{q}$, $\mathbf{\dot{c}}$, $\mathbf{\dot{q}}$, $\mathbf{c}_{tg}$). As
in the case of the GRL approach, the initial state at the start of each
episode is set at the nominal joint pose ($\mathbf{c}_{0},\mathbf{q}_{0}$),
with zero velocity. We also encourage smoothness by penalizing the quantity
$\mathbf{\tilde{q}}^{j}-\mathbf{\tilde{q}}^{j-1}$. Because of the different
units, to have better conditioning in the gradient of the NN, we scale each
state variable against its range. With respect to the GRL approach we increase
the batch size to 512 to collect more observations for training. We provide
the feasibility rewards $C_{i}$ at each loop as differential w.r.t. the
previous loop. This approach is meant to ease learning while converging to the
same solution [38]. At the end of each episode, we strongly penalize the lack
of a touch-down event, and use the same target function (9) to encourage
landing close to the target. The fact that we provide this quantity at every
step enables us to achieve an informative reward [38] also before the
touchdown. Following the curriculum learning idea [39], we gradually increase
the difficulty of the jump enlarging the bounds of the training region (where
the targets are sampled) in accordance with the number of episodes and the
average reward.

###  IV-A Policy performance: the feasibility region

We tested the agent in inference mode, for omni-directional jumps at different
heights, for 726 target positions uniformly spaced on a grid (test region) of
the same shape of the training region. The test region is 20$\%$ bigger than
the training region, in order to demonstrate the generalisation abilities of
the system.

The policy was periodically evaluated on the test region set in order to
assess the evolution of the models stored during the training phase. To
measure the quality of a jump, we used the Relative Percentual Error (RPE),
which we define as the distance between the touch-down and the desired target
point, divided by the jump length. The feasibility region represents the area
where the agent is capable of an accurate landing, i.e. RPE $\leq 10\%$.

####  IV-A1 Performance baseline: Trajectory Optimisation

We compared the approach with the baseline FDDP approach repeating the same
optimisation for all the points in the test region without changing the
costs’weights and limiting the number of iterations to $500$. For optimal
control, the average computation time was $17$ $s$ for back jumps and $7.6$ s
for front jumps while a single evaluation of the NN requires only $0.7$ $ms$.
Fig. 3 (right) shows that a reasonable accuracy is obtained for landing
locations in the front of the robot, while FDDP behaves poorly for locations
in the back of the robot. Computing the mean RPE separately for the back
region and the front region, we obtained an accuracy of 52 $\%$ and 16.5$\%$,
respectively.

####  IV-A2 Performance of the GRL approach

We repeated the simulations on the test region using the GRL approach with the
default NN and with a NN where we halved the number of neurons in each hidden
layer (half NN). Fig. 4 shows that, in both cases, with our GRL approach the
RPE decreases (i.e. accuracy increases) monotonically with the number of
episodes. In the front case, going from approximately 40$\%$ to 16$\%$ for
100k episodes. A satisfactory level (i.e. RPE 20$\%$ for front jumps) is
already achieved after 50k episodes. All the feasibility constraints turn out
to be mostly satisfied after 10k episodes. The figure also shows that the GRL
approach always outperforms the standard optimisation method in terms of jump
accuracy in the case of back jumps, achieving a comparable accuracy for front
jumps.

Figure 4: Plot of the average RPE as a function of the number of training
episodes.

Halving the number of neurons, the two models behave with similar accuracy,
showing that the RL model could ideally be further simplified. From Fig. 3 we
can observe that the feasibility region expands with the number of training
epochs. In the same figure we can see that the test region (black cylinder) is
bigger than the region where we trained the NN (gray cylinder), demonstrating
its extrapolation capabilities. It is important to remark that the GRL
approach is also capable, to some extent, to learn the dynamics and compensate
the tracking inaccuracies of the underlying low-level controller. This is
shown in the accompanying video 333Link to the accompanying video by the
purple ball that represents the ideal landing location (i.e. if the CoM lift-
off velocity associated to each action was perfectly tracked). The purple ball
is different from the target location (blue ball) because of tracking
inaccuracies but the agent learned to provide a lift-off velocity that
compensates for these, managing to reach accurately the desired location (blue
ball). In the same video we show how that the quality of the jumps steadily
improves with the number of training episodes.

####  IV-A3 Performance baseline: end to end RL

We did not achieve satisfactory results with the E2E approach. The reward
function had an erratic behaviour during the training and even after 1M of
training steps there were only a few targets out of $726$ where the algorithm
managed to attain an error below $10$ $\%$. More details are reported in the
accompanying video.

##  V Conclusions

In this work we proposed a guided RL based strategy to perform omni-
directional jumps on uneven terrain with a legged robot. Exploiting some
domain knowledge and taking a few assumptions on the shape of the jump
trajectory, we have shown that, in a few thousand episodes of training, the
agent obtains the ability to jump in a big area maintaining high accuracy
while reaching the boundary of his performances (deriving from the the
physical limitations of the machine). The approach manages also to learn and
proficiently compensate from the tracking inaccuracies of the low-level
controller. The proposed approach is very efficient (it requires a small
number of training episodes to reach a good performance), it achieves a good
generalisation (e.g., by executing jumps in a region 20% larger than the one
used for training), and it outperforms a standard end-to-end RL that resulted
not able to learn the jumping motion. Compared to optimal control, the GRL
approach 1) achieves the same level of performance in front jumps, but is also
able to perform backward jumps (optimal control is not) 2) requires several
orders of magnitude lower computation time. In the future, we plan to extend
the approach to a full quadruped robot, considering not only linear but also
angular motions. Leveraging the angular part we can build a framework that is
able to perform a variety of jumping motions (e.g. twist, somersault, barrel
jumps) on inclined surfaces. We are also seeking ways to improve robustness by
including robot non idealities in the learning phase and to speed up the
training phase by leveraging parallel computation.

## References

  * [1] F. Jenelten, R. Grandia, F. Farshidian, and M. Hutter, “Tamols: Terrain-aware motion optimization for legged systems,” _IEEE Transactions on Robotics_ , vol. 38, no. 6, pp. 3395–3413, 2022.
  * [2] F. Roscia, M. Focchi, A. D. Prete, D. G. Caldwell, and C. Semini, “Reactive landing controller for quadruped robots,” 2023.
  * [3] H.-W. Park, P. M. Wensing, and S. Kim, “High-speed bounding with the mit cheetah 2: Control design and experiments,” _The International Journal of Robotics Research_ , vol. 36, no. 2, pp. 167–192, 2017. [Online]. Available: https://doi.org/10.1177/0278364917694244
  * [4] J. K. Yim, B. R. P. Singh, E. K. Wang, R. Featherstone, and R. S. Fearing, “Precision robotic leaping and landing using stance-phase balance,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 2, pp. 3422–3429, 2020\.
  * [5] C. Nguyen and Q. Nguyen, “Contact-timing and trajectory optimization for 3d jumping on quadruped robots,” in _2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2022, pp. 11 994–11 999.
  * [6] M. Chignoli and S. Kim, “Online trajectory optimization for dynamic aerial motions of a quadruped robot,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021, pp. 7693–7699.
  * [7] G. García, R. Griffin, and J. Pratt, “Time-varying model predictive control for highly dynamic motions of quadrupedal robots,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , 2021, pp. 7344–7349.
  * [8] M. Chignoli, S. Morozov, and S. Kim, “Rapid and reliable quadruped motion planning with omnidirectional jumping,” in _2022 International Conference on Robotics and Automation (ICRA)_ , 2022, pp. 6621–6627.
  * [9] C. Mastalli, W. Merkt, G. Xin, J. Shim, M. Mistry, I. Havoutis, and S. Vijayakumar, “Agile maneuvers in legged robots:a predictive control approach,” _ArXiv_ , 2022.
  * [10] T. P. Lillicrap, J. J. Hunt, A. Pritzel, N. M. O. Heess, T. Erez, Y. Tassa, D. Silver, and D. Wierstra, “Continuous control with deep reinforcement learning,” _CoRR_ , vol. abs/1509.02971, 2015.
  * [11] C. Gehring, S. Coros, M. Hutler, C. Dario Bellicoso, H. Heijnen, R. Diethelm, M. Bloesch, P. Fankhauser, J. Hwangbo, M. Hoepflinger, and R. Siegwart, “Practice Makes Perfect: An Optimization-Based Approach to Controlling Agile Motions for a Quadruped Robot,” _IEEE Robot. Autom. Mag._ , vol. 23, no. 1, pp. 34–43, 2016.
  * [12] J. Hwangbo, J. Lee, A. Dosovitskiy, D. Bellicoso, V. Tsounis, V. Koltun, and M. Hutter, “Learning agile and dynamic motor skills for legged robots,” _Science Robotics_ , vol. 4, no. 26, pp. 1–20, 2019.
  * [13] X. Peng, E. Coumans, T. Zhang, T.-W. Lee, J. Tan, and S. Levine, “Learning agile robotic locomotion skills by imitating animals,” _Robotics: Science and Systems_ , 07 2020.
  * [14] G. Ji, J. Mun, H. Kim, and J. Hwangbo, “Concurrent Training of a Control Policy and a State Estimator for Dynamic and Robust Legged Locomotion,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 2, pp. 4630–4637, 2022\.
  * [15] N. Rudin, D. Hoeller, P. Reist, and M. Hutter, “Learning to walk in minutes using massively parallel deep reinforcement learning,” in _Conference on Robot Learning_. PMLR, 2022, pp. 91–100.
  * [16] P. Fankhauser, M. Hutter, C. Gehring, M. Bloesch, M. A. Hoepflinger, and R. Siegwart, “Reinforcement learning of single legged locomotion,” _IEEE Int. Conf. Intell. Robot. Syst._ , pp. 188–193, 2013.
  * [17] M. Bogdanovic, M. Khadiv, and L. Righetti, “Model-free reinforcement learning for robust locomotion using demonstrations from trajectory optimization,” _Frontiers in Robotics and AI_ , vol. 9, 2022. [Online]. Available: https://www.frontiersin.org/articles/10.3389/frobt.2022.854212
  * [18] G. Bellegarda and Q. Nguyen, “Robust quadruped jumping via deep reinforcement learning,” _arXiv preprint arXiv:2011.07089_ , 2020.
  * [19] G. Grandesso, E. Alboni, G. P. Papini, P. M. Wensing, and A. D. Prete, “CACTO: Continuous Actor-Critic With Trajectory Optimization - Towards Global Optimality,” _IEEE Robotics and Automation Letters_ , vol. 8, no. 6, pp. 3318–3325, 2023.
  * [20] X. B. Peng and M. van de Panne, “Learning locomotion skills using deeprl: Does the choice of action space matter?” in _Proc. ACM SIGGRAPH / Eurographics Symposium on Computer Animation_ , 2017.
  * [21] G. Bellegarda and K. Byl, “Training in task space to speed up and guide reinforcement learning,” in _2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2019, pp. 2693–2699.
  * [22] S. Chen, B. Zhang, M. W. Mueller, A. Rai, and K. Sreenath, “Learning Torque Control for Quadrupedal Locomotion,” _ArXiv_ , 2022. [Online]. Available: http://arxiv.org/abs/2203.05194
  * [23] M. Aractingi, P.-A. Léziart, T. Flayols, J. Perez, T. Silander, and P. Souères, “Controlling the Solo12 Quadruped Robot with Deep Reinforcement Learning,” Aug. 2022, working paper or preprint. [Online]. Available: https://hal.laas.fr/hal-03761331
  * [24] P. Henderson, J. Hu, J. Romoff, E. Brunskill, D. Jurafsky, and J. Pineau, “Towards the systematic reporting of the energy and carbon footprints of machine learning,” _The Journal of Machine Learning Research_ , vol. 21, no. 1, pp. 10 039–10 081, 2020.
  * [25] A. M. Zador, “A critique of pure learning and what artificial neural networks can learn from animal brains,” _Nature communications_ , vol. 10, no. 1, p. 3770, 2019.
  * [26] H. Shen, J. Yosinski, P. Kormushev, D. G. Caldwell, and H. Lipson, “Learning fast quadruped robot gaits with the rl power spline parameterization,” _Cybernetics and Information Technologies_ , vol. 12, no. 3, pp. 66–75, 2013\. [Online]. Available: https://doi.org/10.2478/cait-2012-0022
  * [27] T. Kim and S.-H. Lee, “Quadruped locomotion on non-rigid terrain using reinforcement learning,” _ArXiv_ , vol. abs/2107.02955, 2021. [Online]. Available: https://api.semanticscholar.org/CorpusID:235755354
  * [28] Y. Ji, Z. Li, Y. Sun, X. B. Peng, S. Levine, G. Berseth, and K. Sreenath, “Hierarchical reinforcement learning for precise soccer shooting skills using a quadrupedal robot,” in _2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2022, pp. 1479–1486.
  * [29] S. Fujimoto, H. Hoof, and D. Meger, “Addressing function approximation error in actor-critic methods,” in _International Conference on Machine Learning_ , 2018, pp. 1582–1591.
  * [30] M. Grzes, _Reward shaping in episodic reinforcement learning_. ACM, 2017.
  * [31] M. Focchi, F. Roscia, and C. Semini, “Locosim: an open-source cross-platform robotics framework,” 2023.
  * [32] R. Budhiraja, J. Carpentier, C. Mastalli, and N. Mansard, “Differential Dynamic Programming for Multi-Phase Rigid Contact Dynamics,” in _IEEE International Conference on Humanoid Robots_ , 2018.
  * [33] C. Mastalli, R. Budhiraja, W. Merkt, G. Saurel, B. Hammoud, M. Naveau, J. Carpentier, L. Righetti, S. Vijayakumar, and N. Mansard, “Crocoddyl: An Efficient and Versatile Framework for Multi-Contact Optimal Control,” _Proceedings - IEEE International Conference on Robotics and Automation_ , no. March 2020, pp. 2536–2542, 2020.
  * [34] J. Carpentier, G. Saurel, G. Buondonno, J. Mirabel, F. Lamiraux, O. Stasse, and N. Mansard, “The Pinocchio C++ library – A fast and flexible implementation of rigid body dynamics algorithms and their analytical derivatives,” in _IEEE International Symposium on System Integrations (SII)_ , 2019.
  * [35] S. Gangapurwala, L. Campanaro, and I. Havoutis, “Learning Low-Frequency Motion Control for Robust and Dynamic Robot Locomotion,” _arXiv_ , 2022\. [Online]. Available: http://arxiv.org/abs/2209.14887
  * [36] T. Z. Zhao, V. Kumar, S. Levine, and C. Finn, “Learning Fine-Grained Bimanual Manipulation with Low-Cost Hardware,” _arXiv_ , 2023. [Online]. Available: http://arxiv.org/abs/2304.13705
  * [37] M. Aractingi, P.-A. Léziart, T. Flayols, J. Perez, T. Silander, P. Souères, and P. Soueres, “Controlling the Solo12 Quadruped Robot with Deep Reinforcement Learning,” _Scientific Reports_ , vol. 13, 2023\.
  * [38] S. H. Jeon, S. Heim, C. Khazoom, and S. Kim, “Benchmarking potential based rewards for learning humanoid locomotion,” in _2023 IEEE International Conference on Robotics and Automation (ICRA)_ , 2023, pp. 9204–9210.
  * [39] Y. Bengio, J. Louradour, R. Collobert, and J. Weston, “Curriculum learning,” _ACM International Conference Proceeding Series_ , vol. 382, no. November, 2009.


# Fragments in symmetric configurations with block size 3

Grahame Erskine11footnotemark: 1
<EMAIL_ADDRESS>Open University, Milton Keynes, UK Terry
Griggs11footnotemark: 1
<EMAIL_ADDRESS>Jozef Širáň22footnotemark: 2
<EMAIL_ADDRESS>Open University, Milton Keynes, UK and Slovak
University of Technology, Bratislava, Slovakia

###### Abstract

We begin the study of collections of three blocks which can occur in a
symmetric configuration with block size 3, $v_{3}$. Formulae are derived for
the number of occurrences of these and it is shown that the triangle, i.e.
abf, ace, bcd is a basis. It is also shown that symmetric configurations
without triangles exist if and only if $v=15$ or $v\geq 17$. Such
configurations containing “many” triangles are also discussed and a complete
analysis of the triangle content of those with a cyclic automorphism is given.

††Mathematics subject classification: 05B30††Keywords: configuration; fragment

##  1 Introduction

In this paper we will be concerned with symmetric configurations with block
size 3. First we recall the definitions. A _configuration_ $(v_{r},b_{k})$ is
a finite incidence structure with $v$ points and $b$ blocks, with the property
that there exist positive integers $k$ and $r$ such that

  1. (i)

each block contains exactly $k$ points;

  2. (ii)

each point is contained in exactly $r$ blocks; and

  3. (iii)

any pair of distinct points is contained in at most one block.

If $v=b$ (and hence necessarily $r=k$), the configuration is called
_symmetric_ and is usually denoted by $v_{k}$. We are interested in the case
where $k=3$. The blocks will also be called _triples_. A configuration is said
to be _decomposable_ or _disconnected_ if it is the union of two
configurations on distinct point sets. We are primarily interested in
indecomposable (connected) configurations, and so unless otherwise noted, this
is assumed throughout the paper.

It is natural to associate two graphs with a symmetric configuration $v_{3}$.
The first is the _Levi graph_ or _point-block incidence graph_ , obtained by
considering the $v$ points and $v$ blocks of a configuration as vertices, and
including an edge from a point to every block containing it. It follows that
the Levi graph is a cubic (3-regular) bipartite graph of girth at least six.
The second graph is the _incidence graph_ , obtained by considering only the
points as vertices and joining two points by an edge if and only if they
appear together in some block. Thus the incidence graph is regular of valency
6 and order $v$.

The above definition of configuration is the classical one, going back to the
nineteenth century and the work of mathematicians such as Kantor [8] and
Martinetti [9]. It is the meaning of the term in the section of the Handbook
of Combinatorial Designs by Gropp [6]. However, in the last century the term
came to be used more widely to mean any (small) collection of blocks which may
appear in a combinatorial structure. This is the sense in which the term is
used for example in chapter 13 of Triple Systems by Colbourn & Rosa [1]. For
the present paper this double usage of the term is particularly unfortunate as
we wish to study the occurrence of (20th century) configurations in (19th
century) configurations. Therefore, for this reason we will define any
collection of blocks or triples which may appear in a symmetric configuration
$v_{3}$ as a _fragment_.

We begin by developing formulae for the number of occurrences of one-, two-
and three-block fragments in symmetric configurations $v_{3}$. A one-block
fragment is of course just a single triple and there are $v$ of these by
definition. There are two two-block fragments: $A_{1}$, a pair of disjoint
triples and $A_{2}$, a pair of intersecting triples. Here, and below, we
follow the terminology used in [1]. It is elementary to derive the
corresponding formulae. Trivially $a_{2}=3v$ and so
$a_{1}=v(v-1)/2-3v=v(v-7)/2$. Both of these formulae are called _constant_
meaning that they depend only on $v$ and so do not vary over all
configurations $v_{3}$, irrespective of the structure of the individual
configurations. Other fragments are called _variable_.

There are five three-block fragments. Omitting set brackets and commas for
simplicity these are $B_{1}$, abc, def, ghi (3-partial parallel class or
3-PPC); $B_{2}$, abc, def, dgh (hut); $B_{3}$, abc, ade, afg (3-star);
$B_{4}$, abc, cde, efg (3-path); $B_{5}$, abf, ace, bcd (triangle). These are
illustrated in Figure 1.

& $B_{1}$ $B_{2}$ $B_{3}$ $B_{4}$ $B_{5}$

Figure 1: The possible three-block fragments

Let the number of occurrences of fragment $B_{i}$ be $b_{i}$, $1\leq i\leq 5$.
The number of triangles in a symmetric configuration $v_{3}$ is equal to the
number of 6-cycles in its Levi graph and is therefore variable. Let $b_{5}=t$.
Trivially $b_{3}=v$. To determine $b_{4}$, consider a pair of intersecting
triples and the eight three-block fragments formed by this pair of triples and
the third block through any of the four points other than the point of
intersection. Then $3v\times 8=6t+2b_{4}$ from which $b_{4}=3(4v-t)$. For
$b_{2}$, again consider a pair of intersecting triples and a third block
through any other point. Then $3v(v-5)\times 3=3t+4b_{4}+6b_{3}+3b_{2}$ from
which $b_{2}=3(v(v-11)+t)$. Finally for $b_{1}$, consider a pair of disjoint
triples and a third block through any other point. Then $v(v-7)(v-6)\times
3/2=b_{4}+4b_{2}+9b_{1}$ from which $b_{1}=(v^{3}-21v^{2}+122v-6t)/6$. For
ease of reference we collect these formulae together in a table.

$b_{1}=(v^{3}-21v^{2}+122v-6t)/6$,
---
$b_{2}=3(v(v-11)+t)$,
$b_{3}=v$,
$b_{4}=3(4v-t)$,
$b_{5}=t$.

From the above, $0\leq t\leq 4v$. It is natural to ask whether symmetric
configurations $v_{3}$ with values of $t$ at either end of the spectrum exist
and how these may be constructed. At the lower end, configurations with no
triangles exist for $v=15$ and all $v\geq 17$. This is proved in Section 2. At
the upper end, the unique configuration $7_{3}$, which is the Fano plane,
contains 28 triangles and is the only configuration attaining the upper bound.
Configurations with many triangles are studied in Section 3. Section 4 deals
with cyclic configurations, i.e. those with a cyclic automorphism. Finally in
Section 5 we discuss some further investigations suggested by the work in this
paper. Computer results on the distribution of the values of $t$ for all
configurations $v_{3}$ for $7\leq v\leq 16$ are given in Table 1 in the
Appendix at the end of the paper.

##  2 Configurations with no triangles

As stated above, the number of triangles in a symmetric configuration $v_{3}$
is equal to the number of 6-cycles in its Levi graph. Therefore in order to
construct configurations with no triangles it is both necessary and sufficient
to construct cubic bipartite graphs with girth equal to 8 (or greater than or
equal to 8). We do this recursively, using the following result.

######  Proposition 2.1.

Let $G$ be a cubic bipartite graph of order $2v$, of girth 8 and which also
contains a 10-cycle. Then there exists a cubic bipartite graph $G^{*}$ of
order $2v+10$, also of girth 8 and containing a 10-cycle.

###### Proof.

Consider the graph $G$. Denote the vertices of the 10-cycle by
$v_{0},v_{1},\ldots,v_{9}$. Let the third vertex to which the vertex $v_{i}$
is adjacent be $w_{i},~{}0\leq i\leq 9$. The vertices $w_{i}$ are distinct.
Further the only two such vertices which can be adjacent are
$\\{w_{i},w_{i+5}\\},~{}0\leq i\leq 4$. Delete the five edges
$\\{v_{2i-1},v_{2i}\\},~{}1\leq i\leq 5$, subscripts modulo 10, to form a
graph $G^{\prime}$. Now construct a graph $H$ consisting of a 10-cycle
$u_{0},u_{3},u_{4},u_{7},u_{8},u_{1},u_{2},u_{5},u_{6},u_{9}$ where in
addition each vertex $u_{i}$ is adjacent to a vertex $v_{i}$. Identify the
vertices $v_{i}$ in the two graphs $G^{\prime}$ and $H$ to form the graph
$G^{*}$. Clearly the graph is cubic, bipartite (as represented by the
black/white colouring of the vertices as shown in Figure 2), of order $2v+10$
and contains a 10-cycle. It remains to prove that it has girth 8.

Figure 2: Graphs from the proof of Proposition 2.1

&
$v_{0}$$v_{1}$$v_{2}$$v_{3}$$v_{4}$$v_{5}$$v_{6}$$v_{7}$$v_{8}$$v_{9}$$w_{5}$$w_{0}$$w_{1}$$w_{2}$$w_{9}$$w_{8}$$w_{7}$$w_{6}$$w_{4}$$w_{3}$
$G$ $G^{\prime}$
$u_{3}$$u_{4}$$u_{7}$$u_{8}$$u_{1}$$u_{2}$$u_{5}$$u_{6}$$u_{9}$$u_{0}$$v_{2}$$v_{3}$$v_{4}$$v_{7}$$v_{0}$$v_{9}$$v_{6}$$v_{5}$$v_{1}$$v_{8}$
$H$ Let $C$ be a cycle of $G^{*}$. If $C$ contains only edges of $G^{\prime}$
or only edges of $H$, then $C$ has length at least 8. If not, then $C$ must
contain at least three edges from $H$ including two edges $\\{u_{m},v_{m}\\}$
and $\\{u_{n},v_{n}\\}$, $0\leq m,n\leq 9$. We need to show that it now
contains at least five edges from $G^{\prime}$. If $C$ contains only one edge
from $G^{\prime}$ it would have to be $\\{v_{2i},v_{2i+1}\\}$ for some $i$
such that $0\leq i\leq 4$. But then $C$ would actually contain seven edges
from $H$ and so have length 8. Clearly $C$ cannot contain only two edges from
$G^{\prime}$. If it contains just three edges then these would have to be
$\\{v_{i},w_{i}\\}$, $\\{w_{i},w_{i+5}\\}$ and $\\{w_{i+5},v_{i+5}\\}$ for
some $i$ such that $0\leq i\leq 4$ and then $C$ would contain five edges from
$H$ and so again have length 8. The remaining possibility is that $C$ contains
just four edges from $G^{\prime}$, in which case vertices $v_{m}$ and $v_{n}$
are in the same partition and their distance in $H$ is at least 4. ∎

We can now prove the main result of this section.

######  Theorem 2.2.

There exists a symmetric configuration $v_{3}$ with no triangles if and only
if $v=15$ or $v\geq 17$.

###### Proof.

The smallest cubic graph of girth 8 has order 30 and is Tutte’s 8-cage. Thus
there are no configurations $v_{3}$ with no triangles for $7\leq v\leq 14$.
The graph is bipartite, unique, contains a 10-cycle and is the Levi graph of
the Cremona-Richmond configuration or generalised quadrangle GQ($2,2$). It has
an elegant construction as follows. Let $S$ be a set of cardinality 6. Define
the points of the configuration to be the set of unordered pairs of elements
of $S$ (i.e. subsets of cardinality 2) and the lines to be the partitions of
$S$ into disjoint pairs. An exhaustive computer search shows that there is no
cubic bipartite graph of order 32 and girth 8, a fact confirmed on page 734 of
[10].

From Proposition 2.1 in order to complete the proof of the theorem, it
suffices to exhibit cubic bipartite graphs of girth 8 also containing a
10-cycle of orders 34, 36, 38 and 42. There is a unique graph of order 34 and
girth 8 and three graphs of order 36 and girth 8, again see page 734 of [10].
We have determined that all these four graphs are bipartite and contain a
10-cycle. The corresponding configurations are given below.

012 034 056 178 19a 2bc 2de 37b 39d 48e 4af 58c 5df 6ab 6eg 7fg 9cg

012 034 056 178 19a 2bc 2de 37b 39d 48c 4af 58g 5ae 6ch 6df 7eh 9gh bfg

012 034 056 178 19a 2bc 2de 37b 39d 48c 4fg 58e 59f 6ch 6dg 7fh abg aeh

012 034 056 178 19a 2bc 2de 37b 39d 48e 4af 57g 5df 68c 69h acg bfh egh

Configurations corresponding to bipartite graphs of orders 38 and 42 of girth
8 and containing a 10-cycle are also given below.

012 034 056 178 19a 2bc 2de 37b 39d 48c 4af 57e 5ag 68h 69i bfh cgi dgh efi

012 034 056 178 19a 2bc 2de 37b 39d 48c 4af 57e 5gh 68i 6fj 9gk ahi bfg cjk
dhj eik

This completes the proof of the theorem. ∎

##  3 Configurations with many triangles

For any symmetric configuration $v_{3}$, say $C$, we let $t(C)$ denote the
number of triangles in $C$. For any $v\geq 7$ we let $T(v)$ denote the maximum
value of $t(C)$ for any configuration $v_{3}$.

######  Theorem 3.1.

For any $v\geq 8$, $T(v)\leq 3v$.

###### Proof.

As noted above, the number of triangles in a configuration is precisely the
number of 6-cycles in its Levi graph. We therefore count the largest possible
number of 6-cycles in a cubic bipartite graph of girth 6. To do this it is
most convenient to count the number of 6-cycles through any given edge.

Referring to Figure 3, the number of 6-cycles through the edge $ab$ is equal
to the number of edges between those vertices at distance 2 from $ab$, i.e.
the 8 vertices between the dotted lines. There are a maximum of 8 such edges;
but if all 8 are present then the graph is the Heawood graph, which is the
Levi graph of the unique configuration $7_{3}$.

So suppose that the maximum number of 6-cycles through any vertex is 7, as in
the figure. Then the 14 vertices at maximum distance 2 from $ab$ form a
Heawood graph with a single edge deleted; since the Heawood graph is edge-
transitive we may assume without loss that the missing edge is $uv$ as in the
figure. Then the subgraph $H$ induced by these 14 vertices is joined to the
rest of the graph by edges $ux$ and $vy$. Now $H$ contains exactly 20
6-cycles, since a Heawood graph has 28 and removing the edge $uv$ reduces this
by 8. Thus the average number of 6-cycles through any edge in $H$ is 6.
Moreover, since $u$ and $v$ are at distance 5 in $H$, there are no 6-cycles
through either edge $ux$ or $vy$.

The graph may contain multiple edges lying in 7 6-cycles, so there may be more
than one of these induced subgraphs isomorphic to $H$, joined to the remainder
of the graph by edges which are in no 6-cycle. The remainder of the edges in
the graph not in any of these induced subgraphs must lie in a maximum of 6
6-cycles, and this includes a non-zero number of edges not in any 6-cycle at
all. Thus if there is at least one edge lying in 7 6-cycles, the average
number of 6-cycles through all the edges in the graph is strictly less than 6.
Since there are $3v$ edges in the Levi graph and each 6-cycle goes through 6
edges, it follows that there are strictly fewer than $3v$ 6-cycles in the Levi
graph and hence triangles in the configuration.

The remaining possibility to achieve equality in the bound is that there are
precisely 6 6-cycles through each edge in the Levi graph. This situation is
attained by the unique configuration $8_{3}$. ∎

$a$$b$$u$$v$$x$$y$ Figure 3: An edge $ab$ contained in 7 6-cycles

We note that the proof of Theorem 3.1 implies that equality exists in the
bound if and only if every edge in the Levi graph lies in exactly 6 6-cycles.
Such a graph is called _edge-girth-regular_ and these objects were studied in
[7]. The only known edge-girth-regular graph of degree 3, girth 6 and with
every edge lying in 6 girth cycles is the Möbius-Kantor graph, which is the
Levi graph of the unique configuration $8_{3}$. It is an open question whether
there exist further configurations where equality is achieved in the bound of
Theorem 3.1.

######  Theorem 3.2.

$3\geq\limsup_{v\to\infty}\frac{T(v)}{v}\geq\liminf_{v\to\infty}\frac{T(v)}{v}\geq\frac{20}{7}.$

###### Proof.

The upper bound for the $\limsup$ follows immediately from Theorem 3.1. To
demonstrate the lower bound, we construct an infinite family of configurations
on $7n$ points, $n\geq 2$, with exactly $20n$ triangles. To do this we take
$n$ copies $H_{1},H_{2},\ldots,H_{n}$ of the edge-deleted Heawood graph $H$
from the proof of Theorem 3.1 and connect them via their joining edges in a
cyclic manner. The construction is illustrated in Figure 4 for the case $n=3$.
Since the joining edges do not lie in any 6-cycles, the resulting graph is the
Levi graph of a configuration with $7n$ points and $20n$ triangles as
required. To complete the proof we need to show that for any $\epsilon>0$,
there is some $v_{0}$ such that for any $v\geq v_{0}$ we can construct a
configuration $C$ on $v$ points with $t(C)/v\geq 20/7-\epsilon$. It is easily
seen that we can do this by constructing a graph with a large number $n$ of
copies of $H$ plus one edge-deleted Levi graph of a configuration on
$7,8,9,10,11,12$ or $13$ points, and join them cyclically as before. We can
choose $n$ large enough so that the number of triangles in the resulting
configuration $v_{3}$ will be arbitrarily close to $20v/7$. ∎

$H_{1}$$H_{2}$$H_{3}$ Figure 4: The Levi graph of a configuration on 21 points
with 60 triangles

We note that the graph depicted in Figure 4 has appeared before in the
literature. For example, in previous works of the present authors it arises in
connection with colouring problems [2] and embeddings [3].

##  4 Cyclic configurations

A symmetric configuration $v_{3}$ is cyclic if it admits an automorphism of
order $v$. Such a configuration can be realised on the set
$\\{0,1,\ldots,v-1\\}$ as a single orbit under the mapping $i\mapsto i+1$ (mod
$v$) of a starter block
$\\{0,a,a+b\\},~{}a,b\in\\{1,2,\ldots,v-2\\},~{}a+b\leq v-1$. The orbit can be
described as a cyclically ordered triple $\langle a,b,c\rangle$ where
$a+b+c=v$. However in order to generate a configuration, certain requirements
need to be met: $1\leq a,b,c\leq v-3,~{}a\neq b\neq c\neq a$ and if $v$ is
even, $a,b,c\neq v/2$. Further, for the configuration to be connected we must
have $\gcd(a,b,c)=1$. We are able to give a complete analysis of the
occurrence of triangles in these configurations.

######  Theorem 4.1.

Let $\mathcal{C}$ be a cyclic symmetric configuration $v_{3}$. Then the number
of triangles contained in $\mathcal{C}$ is (i) $4v$ if $v=7$, (ii) $3v$ if
$v=8$, (iii) $7v/3$ if $v=9$ and (iv) $v,4v/3$ or $2v$ if $v\geq 10$.

###### Proof.

The configurations $7_{3}$ and $8_{3}$ are unique up to isomorphism and are
cyclic. It is easily determined that the number of triangles they contain are
28 and 24 respectively. Of the three non-isomorphic configurations $9_{3}$,
only one is cyclic and again it is easily determined that the number of
triangles is 21.

Now let $v\geq 10$. The incidence graph is a Cayley graph
Cay($\mathbb{Z}_{v},S$) where the connection set $S=\\{\pm a,\pm b,\pm c\\}$.
Triangles in the configuration will occur as 3-cycles in the incidence graph.
We count these by identifying edges which sum to zero, noting that $(x,y,z)$,
$(y,z,x)$, $(z,x,y)$, $(-x,-z,-y)$, $(-z,-y,-x)$ and $(-y,-x,-z)$ all count
the same triangles. First consider the two sets of edges $(a,b,c)$ and
$(a,c,b)$; the first of these correspond to blocks of the configuration and
the second to $v$ triangles in the configuration.

There are two further situations where triangles occur. The first of these is
where, without loss of generality, either $b=2a$ or $c=2a$, i.e. the
cyclically ordered triple $\langle a,b,c\rangle$ is either $\langle
a,2a,v-3a\rangle$. or $\langle a,v-3a,2a\rangle$. (Note that $b=v-2a$ and
$c=v-2a$ are not possible.) In this case there are a further $v$ triangles
corresponding to edges $(a,a,-2a)$ in the incidence graph. Potentially an
orbit of this form may generate further triangles from 3-cycles as follows;
$(2a,2a,3a)$, $(3a,3a,a)$ and $(3a,3a,2a)$. If $\gcd(a,v)=1$, the first two
can only appear when $v=7$ and the third when $v=8$. This accounts for the
“extra” triangles in these configurations noted above. If $\gcd(a,v)>1$, the
configurations are disconnected.

The second situation can only occur when $v\equiv 0$ (mod 3) and is without
loss of generality when $a=v/3$ or $a=2v/3$. There are $v/3$ triangles
$\langle v/3,v/3,v/3\rangle$ in the configuration corresponding to edges
$(v/3,v/3,v/3)$ in the incidence graph. It remains to consider whether both of
these two situations can occur simultaneously. There are six possibilities:
$(v/3,v/6,v/2)$, $(v/3,v/2,v/6)$, $(v/3,2v/9,4v/9)$, $(v/3,4v/9,2v/9)$,
$(2v/3,v/9,2v/9)$ and $(2v/3,2v/9,v/9)$. The first two can only occur when
${v=6s,s>1}$ and the configurations will be disconnected. The other four occur
when ${v=9s,s\geq 1}$. When $s=1$, again this accounts for the “extra”
triangles in the cyclic $9_{3}$ configuration. When $s>1$, the configurations
are disconnected. ∎

##  5 Conclusion

In this paper we have begun the study of so-called fragments in symmetric
configurations with block size 3. This follows from similar work with respect
to Steiner triple systems which was published nearly 30 years ago [4]. An
entire chapter of [1] is devoted to this and many other aspects of this work
due to various authors. We are of the opinion that there is still much which
can be done in this area relating to configurations and below we indicate some
possible avenues for future research.

As we have shown, the smallest fragment which can be avoided in an infinite
number of symmetric configurations is the triangle. Such a configuration on 15
points, the Cremona-Richmond configuration, was already known and we have
shown that there exist such configurations $v_{3}$ for all $v\geq 17$. We have
also considered symmetric configurations containing “many” triangles. A
further aspect of this work is that of decomposition. It is relevant to ask
about symmetric configurations whose blocks can be decomposed into triangles
(with one or two extra blocks if $v$ is not divisible by 3). Indeed the whole
question of the decomposition of configurations $v_{3}$ into _any_ given
fragment appears to be open.

The work in this paper may also be extended to non-symmetric configurations
$(v_{r},b_{3})$. Elementary calculations give $b=rv/3$, $a_{2}=r(r-1)v/2$ and
$a_{1}=rv(rv-9r+6)$ for the numbers of blocks, intersecting triples and
disjoint triples respectively; thus an additional parameter $r$, the
replication number, is introduced into the equations. Even further,
configurations $(v_{r},b_{3})_{\lambda}$ where any pair of distinct points is
contained in at most $\lambda$ blocks, can be considered. Those configurations
with $\lambda=2$ are called _spatial_. Symmetric spatial configurations are
the subject of a paper by Gropp [5] and even though this class of
configurations is quite restrictive, results would still be of interest. Two
further two-block fragments, $A_{3}$, abc, abd and $A_{4}$, abc, abc (repeated
triple) may also occur; so, unlike the case where $\lambda=1$, two-block
fragments would not be constant.

However perhaps the most obvious investigation would be to extend the study of
fragments in symmetric configurations $v_{3}$ to those with four blocks. There
are 16 four-block fragments which are illustrated on page 210 of [1] together
with their standard labelling. One of these $C_{7}$, abc, ade, afg, ahi
(4-star) cannot occur. Nevertheless a complete analysis of the number of
occurrences of the remaining 15 fragments would seem to be a lengthy and
perhaps tedious although interesting undertaking. But a subproject could be to
consider only those symmetric configurations which have no triangles, i.e.
those whose Levi graph has girth at least 8. This would eliminate $C_{6}$,
$C_{8}$, $C_{11}$, $C_{12}$, $C_{14}$, $C_{15}$ and $C_{16}$ from the
configuration leaving just 8 fragments to consider which is a much more
feasible investigation. Within this scenario an 8-cycle in the Levi graph
corresponds to a square, i.e. fragment $C_{10}$, abx, bcy, cdz, daw. Thus the
construction of cubic bipartite graphs of girth greater than or equal to 10
would give configurations $v_{3}$ avoiding $C_{10}$ and all fragments
containing triangles. A result analogous to Proposition 2.1 would go a long
way to establishing this. There is still much to investigate and we hope that
the reader will be inspired to consider some of these ideas.

## Appendix

$v$ | $t$ | Count
---|---|---
7 | 28 | 1
8 | 24 | 1
9 | 18 | 1
9 | 20 | 1
9 | 21 | 1
10 | 17 | 2
10 | 18 | 3
10 | 19 | 2
10 | 20 | 3
11 | 15 | 1
11 | 16 | 10
11 | 17 | 7
11 | 18 | 7
11 | 19 | 3
11 | 20 | 1
11 | 21 | 1
11 | 22 | 1
12 | 12 | 1
12 | 13 | 8
12 | 14 | 22
12 | 15 | 48
12 | 16 | 60
12 | 17 | 41
12 | 18 | 24
12 | 19 | 14
12 | 20 | 5
12 | 21 | 3
12 | 22 | 1
12 | 24 | 2
13 | 9 | 1
13 | 10 | 2
13 | 11 | 12
13 | 12 | 67
13 | 13 | 190
13 | 14 | 371
13 | 15 | 418

$v$ | $t$ | Count
---|---|---
13 | 16 | 409
13 | 17 | 265
13 | 18 | 156
13 | 19 | 74
13 | 20 | 37
13 | 21 | 14
13 | 22 | 9
13 | 23 | 4
13 | 24 | 3
13 | 25 | 1
13 | 26 | 1
13 | 28 | 1
13 | 32 | 1
14 | 6 | 1
14 | 8 | 5
14 | 9 | 24
14 | 10 | 145
14 | 11 | 521
14 | 12 | 1512
14 | 13 | 2901
14 | 14 | 4086
14 | 15 | 4121
14 | 16 | 3247
14 | 17 | 2236
14 | 18 | 1304
14 | 19 | 640
14 | 20 | 335
14 | 21 | 159
14 | 22 | 69
14 | 23 | 33
14 | 24 | 36
14 | 25 | 5
14 | 26 | 2
14 | 27 | 4
14 | 28 | 6
14 | 31 | 2

$v$ | $t$ | Count
---|---|---
14 | 32 | 3
14 | 40 | 1
15 | 0 | 1
15 | 4 | 3
15 | 6 | 12
15 | 7 | 40
15 | 8 | 254
15 | 9 | 1129
15 | 10 | 4252
15 | 11 | 11877
15 | 12 | 24510
15 | 13 | 38017
15 | 14 | 44834
15 | 15 | 41585
15 | 16 | 32177
15 | 17 | 20914
15 | 18 | 12585
15 | 19 | 6517
15 | 20 | 3341
15 | 21 | 1617
15 | 22 | 851
15 | 23 | 404
15 | 24 | 206
15 | 25 | 62
15 | 26 | 54
15 | 27 | 39
15 | 28 | 23
15 | 29 | 6
15 | 30 | 12
15 | 31 | 9
15 | 32 | 6
15 | 33 | 1
15 | 38 | 2
15 | 40 | 1
16 | 3 | 1
16 | 4 | 5

$v$ | $t$ | Count
---|---|---
16 | 5 | 40
16 | 6 | 277
16 | 7 | 1699
16 | 8 | 8325
16 | 9 | 31782
16 | 10 | 92432
16 | 11 | 206506
16 | 12 | 357339
16 | 13 | 480580
16 | 14 | 517343
16 | 15 | 458294
16 | 16 | 344331
16 | 17 | 225866
16 | 18 | 133408
16 | 19 | 72369
16 | 20 | 37472
16 | 21 | 18591
16 | 22 | 9255
16 | 23 | 4570
16 | 24 | 2123
16 | 25 | 871
16 | 26 | 650
16 | 27 | 310
16 | 28 | 169
16 | 29 | 101
16 | 30 | 82
16 | 31 | 43
16 | 32 | 18
16 | 33 | 5
16 | 34 | 1
16 | 35 | 1
16 | 36 | 12
16 | 37 | 4
16 | 38 | 1
16 | 39 | 1

Table 1: Numbers of triangles $t$ in configurations $v_{3}$

## Acknowledgements

The third author acknowledges support from the APVV Research Grants 15-0220
and 17-0428, and the VEGA Research Grants 1/0142/17 and 1/0238/19.

## References

  * [1] C. J. Colbourn and A. Rosa.  Triple Systems.  Clarendon Press, New York (1999).
  * [2] G. Erskine, T. S. Griggs and J. Širáň.  Colouring problems for symmetric configurations with block size 3.  J. Combin. Des., 29(6) (2021), 397–423.
  * [3] G. Erskine, T. S. Griggs and J. Širáň.  On the upper embedding of symmetric configurations with block size 3.  Discrete Math., 343(4) (2020), article no. 111774.
  * [4] M. J. Grannell, T. S. Griggs and E. Mendelsohn.  A small basis for four-line configurations in Steiner triple systems.  J. Combin. Des., 3 (1994), 51–59.
  * [5] H. Gropp.  On symmetric spatial configurations.  Discrete Math., 125 (1994), 201–209.
  * [6] H. Gropp.  Configurations.  Handbook of Combinatorial Designs 2/e, C. J. Colbourn and J. H. Dinitz (eds.), CRC Press, Boca Raton (2007), 353–355.
  * [7] R. Jajcay, G. Kiss and Š. Miklavič.  Edge-girth-regular graphs.  Europ. J. Combin. 72 (2018), 70–82.
  * [8] S. Kantor.  Die Configurationen $(3,3)_{10}$.  Sitzungsber. Akad. Wiss. Wien, Math.-Nat.Kl. 84 (1881), 1291–1314.
  * [9] V. Martinetti.  Sulle configurazioni piane $\mu_{3}$.  Ann. Mat. Pura Appl. 15 (1887), 1–26.
  * [10] G. F. Royle.  Graphs and Multigraphs.  Handbook of Combinatorial Designs 2/e, C. J. Colbourn and J. H. Dinitz (eds.), CRC Press, Boca Raton (2007), 731–740.