MODELING, SIMULATION, AND OPTIMIZATION OF A MONOPOD HOPPING ON YIELDING
TERRAIN

By

Juntao He

Master’s thesis

Submitted in Partial Fulfillment of the

Requirements for the Degree of

MASTER OF SCIENCE IN MECHANICAL ENGINEERING

EVANSTON, Illinios

August 2021

ACKNOWLEDGEMENTS

First of all, I would like to thank Prof. Kevin Lynch, Prof. Paul Umbanhowar
and Dan Lynch for their invaluable help and support on my master’s thesis.

I would also like to thank all the staff at Northwestern University,
especially those who fight on the Covid-19 front line. Without them, I could
never finish my thesis and course of study in these tough times.

Finally, I would like to thank my parents for their support of my master’s
study. They are my strong backing in the Covid pandemic.

###### Table of Contents

  1. 1 Introduction and Background
    1. 1.1 Chrono simulation environment
    2. 1.2 Thesis outline
  2. 2 Chrono environment validation
    1. 2.1 Simulation setup
    2. 2.2 Comparison between DEM results and experiments
    3. 2.3 The influence of the intrusion speed on ground reaction force
  3. 3 Ground Reaction Force model
    1. 3.1 Simulation setup
    2. 3.2 Resistive forces fitting functions
  4. 4 Periodic Hopping of a Monopod
    1. 4.1 Soft ground model
    2. 4.2 Robot model
    3. 4.3 Robot control
      1. 4.3.1 Feedforwrd control
      2. 4.3.2 Feedforward plus feedback control
    4. 4.4 Multiple hops
  5. 5 Conclusion
  6. Bibliography

###### List of Tables

  1. 2.1 Chrono simulation parameters
  2. 2.2 DEM stress gradient results of three different penetration velocities compared to experiments
  3. 3.1 Chrono simulation parameters
  4. 3.2 Coefficients for 2nd order Fourier filtered force functions.
  5. 4.1 Physical properties of the monopod robot.
  6. 4.2 Boundary conditions

###### List of Figures

  1. 2.1 Sketch of a rigid plate intruder
  2. 2.2 Stress results of three different penetration velocities
  3. 2.3 Comparison of experimental and DEM results
  4. 2.4 Normalized error map
  5. 3.1 Bed cooling graph
  6. 3.2 Logarithmic bed cooling graph
  7. 3.3 Raw stress gradient colormaps and their Fourier filtered colormaps
  8. 3.4 1st order $\alpha_{z,x}$ error maps
  9. 3.5 2nd order $\alpha_{z,x}$ error maps
  10. 3.6 3rd order $\alpha_{z,x}$ error maps
  11. 3.7 Stress gradient error bars
  12. 4.1 Physical models of the monopod robot and soft ground.
  13. 4.2 Sketch of phase transitions
  14. 4.3 Planned single-hop monopod trajectories and their control signals
  15. 4.4 Block diagram of feedforward plus feedback control for stance phase
  16. 4.5 CoM trajectories from optimal(desired) and for the two different control strategies
  17. 4.6 Position and velocity errors of two different control strategies for multiple impact velocities
  18. 4.7 Experimental errors vs. foot size
  19. 4.8 Center of mass trajectories for five hops
  20. 4.9 Error of robot CoM trajectories after five hops for different impact velocities
  21. 4.10 Screenshot of demo video

## CHAPTER 1Introduction and Background

Mobile robots have great potential to improve search and rescue, disaster
response, factory inspection and package delivery. Legged locomotion has
significant advantages over wheeled locomotion on yielding surfaces, such as
soil, sand, snow and gravel. However, researchers have mainly focussed on
generating stable walking and running gaits on rigid flat [1] or uneven
terrain [2]. There are far fewer studies about legged locomotion on yielding
ground.
Researchers already made robots to accomplish some locomotion tasks on soft
ground. Earlier work used fixed kinematics (fixed gaits) to move on soft
ground [3] [4]. Only recently has the community started to look at feedforward
control [5]. In [5], Hubicki et al. present a model and fast optimization
formulation which generates accurate motion plans on granular media with
relatively short solving times. Reference [6] achieved flat-footed bipedal
walking on deformable granular terrain with feedback control, but they did so
by seeking to avoid ground penetration. Reference [7] suggests an active
damping controller to reduce the energy cost for monopedal vertical hopping on
soft ground, while reference [8] seeks to minimize energy loss on soft-landing
using feedforward force control as well as optimal impedance.
In this work, I use feedforward control to develop a periodic gait on soft
ground using a ground reaction force model. I test this approach by performing
discrete element method (DEM) simulations of the monoped robot hopping on
granular materials. To stabilize the target gait, I introduce a feedback term.

### 1.1 Chrono simulation environment

To test the force control approach I develop for periodic hopping on soft
ground, experimental validation using a real robot on real ground would be
preferable. However, due to the challenges of working in the lab during the
Covid-19 pandemic, I instead chose to validate the control approach in a
virtual environment, namely Chrono.
Chrono111http://api.projectchrono.org/5.0.0/ is a physics engine with
specialized support for simulating interaction between rigid bodies and
granular media. Chrono supports C++ and Python, but the support for Python is
still under construction, so to take full advantage of Chrono’s features, C++
is recommended. Chrono’s discrete element method (DEM) code works with a GPU
but there is no need for users to write GPU code because Chrono’s developers
already provided libraries for the GPU part. Chrono provides interfaces with
CAD programs like SOLIDWORKS to simplify model building in Chrono. Chrono also
has better integration of rigid-body dynamics with DEM simulation than other
DEM simulators such as LAMMPS222https://lammps.sandia.gov/ and
LIGGGHTS333https://www.cfdem.com/liggghts-open-source-discrete-element-method-
particle-simulation-code.
In the following chapters, I use Chrono to generate resistive forces
approximation functions of the ground reaction force and test the analytic
feedforward force control solution generated by MATLAB using the RFT model.
Advantages of the Chrono simulation environment over real world experiments
include:

  1. 1.

To the extent that the DEM simulation is physically accurate, in-simulation
validation provides more feedback than hardware validation, since potentially
any variable of interest can be measured and recorded. This results in a
better-informed motion planning framework.

  2. 2.

The cost of experimental setup and the cost of experimental failure are both
lower than with hardware experiments.

  3. 3.

The performance of the robot and its control on different soft substrates can
be investigated by simply changing the parameters related to the physical
properties of the granular materials.

### 1.2 Thesis outline

The structure of the thesis is as follows:
Chapter 2 validates the reliability and accuracy of the Chrono simulation
environment by comparing simulation results with experimental results from Li
et al. [3]. Chapter 3 generates functions to predict ground reaction forces
from granular media using a Fourier representation matched to Chrono
simulation data. Chapter 4 develops a control strategy which enables the 1D
hopping monopod to jump back to its original maximum height. Chapter 5
compares the performance of feedforward controls in the DEM simulation for a
range of hop heights and for hoppers with different foot areas. Chapter 6
presents the conclusions.

## CHAPTER 2Chrono environment validation

To validate the Chrono simulation environment, I compared my simulation
results with the experimental results in Li et al. [3].
Refrence [3] finds that the granular resistive stresses are linear with depth,
but only for low intrusion speeds (less than 0.5 m/s), only in the bulk, i.e.,
away from the container sides (to avoid Janssen effects [9]) and at depths
both large enough to avoid nonlinearities due to free particles at the surface
and small enough to avoid nonlinearities due to compression against the
container bottom. For a rigid plate intruder moving in granular media [3], the
lithostatic stresses can be modeled as

$\sigma_{z,x}(|z|,\beta,\gamma)=\left\\{\begin{array}[]{rl}\alpha_{z,x}(\beta,\gamma)|z|&z<0,\\\
0&z>0,\end{array}\right.$ (2.1)

where $\sigma_{z,x}$ are the vertical ($z$-component) and horizontal
($x$-component) resistive stresses, $z$ is the depth of the intruder, and
$\alpha_{z,x}$ are vertical and horizontal stresses per unit depth in the
linear response regime. Note that when the intruder is above the granular
domain ($z>0$), there are no resistive stresses. Angle of attack $\beta$ and
angle of intrusion $\gamma$ are shown in Figure 2.1.

### 2.1 Simulation setup

Figure 2.1 shows a planar view of the plate penetrating into the granular
domain. The black rectangle represents the plate while the colored area
represents the granular media. The plate has dimensions of 3.81 cm $\times$
2.54 cm $\times$ 0.64 cm (area = 9.68 $\textrm{cm}^{2}$) and a density of 2.7
$\mbox{g/cm}^{3}$, which match the dimensions and density of the flat plate
used by Li et al. in [3]. All the spheres in the granular domain have
identical physical properties, specifically a radius of 3.2 mm and a density
of 2.6 $\mbox{g/cm}^{3}$. The spheres are contained in a box with cross
section 24 cm $\times$ 22 cm, and fill the box to a depth of approximately 18
cm. The bed is prepared by filling the container from the bottom to the top
layer by layer. Li et al. used a few different granular materials in their
experiments. I chose 3.2 mm diameter glass spheres instead of those materials
with small radius here because I can use far fewer granular spheres for my
simulation than smaller radius. This makes our simulation time cost more
acceptable. That is, the simulation time for a specific $\gamma-\beta$ pair
intrusion trial is 2 to 3 minutes. For materials with larger radii, fewer
particles contact a fixed size plate, resulting in larger force fluctuations
[10], a topic I return to in the context of feedback control in Chapter 4.
Considering both reducing force fluctuations and simulation time, 3.2 mm is a
good choice for the granular particle radius. Table 2.1 lists the simulation
parameters, for more details, refer to the JSON file in my GitHub repository
111https://github.com/HappyLamb123/Foot-
ROBOT/blob/master/demo_code/plate/demo_GRAN_plate.json. The friction is chosen
the same as used in [3]. Damping coefficients are not mentioned in [3], I set
them as reasonable values for glass.

Figure 2.1:  The rigid plate intrudes into the granular media with velocity $V$, attack angle $\beta$ and intrusion angle $\gamma$. The distance between the CoM of the plate and the granular bed top, $\mid z\mid$, is the intrusion depth.  Time step | 1e-5 s
---|---
Plate dimension | 3.81 cm $\times$ 2.54 cm $\times$ 0.64 cm
Plate density | 2.7 $\mbox{g/cm}^{3}$
Sphere radius | 3.2 mm
Sphere density | 2.6 $\mbox{g/cm}^{3}$
Domain size | 24 cm $\times$ 22 cm $\times$ 18 cm (depth)
Sphere-sphere normal contact stiffness | 1e8 N/m
Sphere-wall normal contact stiffness | 1e8 N/m
Sphere-sphere normal damping coefficient | 6000 N s/m
Sphere-wall normal damping coefficient | 6000 N s/m
Sphere-sphere tangent contact stiffness | 1e6 N/m
Sphere-wall tangent contact stiffness | 1e6 N/m
Sphere-sphere tangent damping coefficient | 0
Sphere-wall tangent damping coefficient | 0
Sphere-sphere static friction coefficient | 0.63
Sphere-wall static friction coefficient | 0.63
Table 2.1: Parameters for Chrono granular simulations.

To reproduce the stress parameters results for 3 mm glass spheres from the
experiments shown in Fig. S4(I, J) in the supplementary materials of [3], I
did 271 simulations with attack angle $\beta$ and intrusion angle $\gamma$
varied from -90 degrees to 90 degrees (three repeated trials for each $\beta$
and $\gamma$ pair). For the intrusion speed $V$, I chose 3 cm/s because it
guarantees both accuracy(how close the simulation results are to experimental
results in [3]) and time-efficiency. Further discussion of the choice of
intrusion speed can be found at the end of this chapter. In each simulation,
$\beta$ and $\gamma$ are fixed during penetration and the absolute value of
intrusion speed $V$ is 3 cm/s. In the Chrono simulation, I record the
horizontal and vertical resistive forces of the granular media every
microsecond into a CSV file. Then, I divide the resistive forces by the area
of the rigid plate to get the stress in the $x$\- and $z$\- direction. To get
the stress gradient, I linearly fit the horizontal and vertical direction
stress from 2 cm to 7.5 cm. The resulting slopes of the lines are the stress
gradient, $\alpha_{z,x}$. I choose the depth range for fitting in order to
avoid surface effects and effects from the container bottom.

### 2.2 Comparison between DEM results and experiments

As [3] suggests, the stresses in the horizontal and vertical directions are
proportional to the penetration depth as Figure 2.2 shows. By linearly fitting
the stress versus depth curves in MATLAB, the stress gradients ($\alpha_{x}$
and $\alpha_{z}$) can be calculated as the slopes of those fitted curves.
Figure 2.3 shows the Chrono simulation results and experimental results of $x$
and $z$ direction stress gradients for different $\beta$ and $\gamma$. These
DEM simulation results show that horizontal and vertical resistive stresses
scale linearly with intruder depth, indicating that the DEM simulations are a
viable alternative to experiments performed with real granular media. The
upper two plots in Figure 2.3 show experimental results from Li et al., while
the lower two are DEM results generated by Chrono. As Figure 2.3 indicates,
the simulation results closely match the experimental results. To better
characterize the accuracy of the Chrono simulation environment, a normalized
error map is presented in Figure 2.4. The normalized errors in Figure 2.4. are
calculated as,

$\varepsilon_{i}=\frac{\alpha_{i,\mathrm{DEM}}-\alpha_{i,\mathrm{Exp}}}{\bar{\alpha}_{i,\mathrm{Exp}}},$
(2.2)

where $i$ is $x$ or $z$, $\alpha_{i,\mathrm{DEM}}$ represents the horizontal
($x$) or vertical ($z$) direction stress per cubic centimeter of DEM results
while $\alpha_{i,\mathrm{Exp}}$ represents experimental results of Li et
al.[3]. The $\bar{\alpha}_{i,\mathrm{Exp}}$ term in the denominator represents
the mean value of all the experimental results .
As the color map in Figure 2.3 indicates, the normalized errors are close to
zero in most regions of the $\beta-\gamma$ plane. Only a few parts are dark
red or blue where the simulation results and experimental results show larger
differences. For $z$-direction normalized errors, the darkest color can be
found in the column where $\gamma$ = 90 degrees. As for $x$-direction
normalized errors, the darkest color can be found in the columns where
$\gamma$ equals 18.435, 45 and 55.31 degrees. Other regions of the
$\beta-\gamma$ plane have very small error values (less than 0.2). So the
simulation results have a good match to the experimental results.

Figure 2.2: Vertical ($z$) and horizontal ($x$) stresses vs. penetration depth
and corresponding linear fits (lines). Similar results for intrusion
velocities of 1, 3, and 5 cm/s indicate that 3 cm/s is a good choice for
balancing simulation time and accuracy relative to experimental results in
[3].  Figure 2.3: Experimental and DEM results for vertical and horizontal
stresses indicate good quantitative agreement for all $\gamma$ and $\beta$.
Figure 2.4: Normalized error map (see equation 2.2).

### 2.3 The influence of the intrusion speed on ground reaction force

In my simulations, the intrusion speed is 3 cm/s while Li et al. used 1 cm/s
in their experiments. I opted for a faster intrusion speed because the rigid
plate can reach the same depth in less time with a higher intrusion speed. A
more acceptable time cost can help us to get more datasets in the same time
period. For the robot control discussion in Chapter 4, a more efficient
simulation time can help reduce the barrier to testing new control policies.
With less time spent in simulation, I can iterate faster and converge to an
effective control policy sooner. To determine the influence of speed
variations, I ran three simulations with different penetration velocities
(Figure 2.2). In those simulations, the simulation parameters are the same as
I discussed above and the rigid plate moves downward in the vertical direction
with a constant speed. The orientation of the plate was also fixed
($\beta$=0). Figure 2.2 plots $\sigma_{x}$ and $\sigma_{z}$ versus penetration
depth $z$ for three constant intrusion velocities $V$. Lines are fitted over a
depth range of $2\ \mathrm{cm}\leq z\leq 7.5\ \mathrm{cm}$. Table 3.1
illustrates the actual mean values of three repeated trials of $\alpha_{x}$
and $\alpha_{z}$ for three velocities and the experimental results in Li’s
paper [3] are added to the last row for comparison. As the table shows, 1 cm/s
is within three percent of Li’s results [3] for $\alpha_{z}$, but increasing
the speed from 1 cm/s to 5 cm/s only changes $\alpha_{z}$ by seven percent.
So, considering the time cost and the accuracy of my simulation, I choose 3
cm/s as the penetration speed for the DEM simulations.

Impact velocity |  $\alpha_{x}$ $[\mathrm{N}/\mathrm{cm}^{3}]$ |  $\alpha_{z}$ $[\mathrm{N}/\mathrm{cm}^{3}]$
---|---|---
Chrono (DEM) results |  |
1 cm/s | 0.004631 | 0.3022
3 cm/s | 0.004155 | 0.2851
5 cm/s | -0.01086 | 0.3232
Experimental results |  |
1 cm/s | 0.0025033 | 0.2931
Table 2.2: Vertically-constrained ($\beta=0,\gamma=\pi/2$) penetration results
for three different velocities. Variation of $\alpha_{x}$ and $\alpha_{z}$ in
three repeated trials is less than 0.003 $\mathrm{N}/\mathrm{cm}^{3}$ in all
cases.

## CHAPTER 3Ground Reaction Force model

Modelling the response of yielding terrains to foot contact is challenging.
The resistive forces from the ground are dependent on the granular
compaction[11] [12], the intruder kinematics (penetration depth and speed) [3]
[13] and the intruder morphology. Reference [13] proposes that the resistive
forces consist of a hydrodynamic-like term and a hydrostatic-like term, and
investigates rapid intrusion by objects that change shape (self-deform)
through passive and active means. [3] presents a granular resistive force
theory to predict the forces on objects intruding relatively slowly(where
inertial effects are negligible) with different directions and orientations.
As discussed in the second paragraph of Chapter 2, [3] proposes that the
granular resistive stresses are linear with depth. To obtain their stress
results, Li et al.[3] intruded a rigid thin plate into a granular bed with
different plate orientaions and movement directions. And then, they got stress
gradients $\alpha_{z,x}(\beta,\gamma)$ (see Eq.2.1) by linearly fitting stress
results at three different depth (2.54 cm, 5.08 cm, and 7.62 cm). Finally,
they performed a discrete Fourier transform on $\alpha_{z,x}(\beta,\gamma)$
results over $-\pi/2<\beta<\pi/2$ and $-\pi<\gamma<\pi$ to obtain a fitting
function. In this chapter, I do something similar but use more modes to get a
more accurate representation. The fitting results in this Chapter will be used
for further locomoting hopper simulations of our group.

### 3.1 Simulation setup

The simulation setup in this chapter is identical to in Figure 2.1, which
shows the planar view of the plate penetrating into the granular domain. The
rigid plate is 5 cm $\times$ 5 cm $\times$ 0.5 cm (area = 25
$\textrm{cm}^{2}$) and the plate mass is 0.25 kg. More details about the
simulation parameters can be found in Table 3.1 and the JSON file in my GitHub
repository 111https://github.com/HappyLamb123/Foot-
ROBOT/blob/master/demo_code/generate_formula/intrude/demo_GRAN_plate.json.

Time step | 5e-5s
---|---
Plate dimension | 5 cm $\times$ 5 cm $\times$ 0.5 cm
Plate density | 20 $\mbox{g/cm}^{3}$
Sphere radius | 3 mm
Sphere density | 2.6 $\mbox{g/cm}^{3}$
Domain size | 40 cm $\times$ 30 cm $\times$ 30 cm
Sphere-sphere normal contact stiffness | 1e8 N/m
Sphere-wall normal contact stiffness | 1e8 N/m
Sphere-sphere normal damping coefficient | 500 N s/m
Sphere-wall normal damping coefficient | 500 N s/m
Sphere-sphere tangent contact stiffness | 1e8 N/m
Sphere-wall tangent contact stiffness | 1e8 N/m
Sphere-sphere tangent damping coefficient | 500
Sphere-wall tangent damping coefficient | 500
Sphere-sphere static friction coefficient | 0.385
Sphere-wall static friction coefficient | 0.385
Table 3.1: Parameters for Chrono granular simulations in Chapter 3.

In Chapter 2, the granular particle container is filled by particles layer by
layer from bottom to top, and no further action is taken for the bed
preparation. A slightly different filling procedure is used in this chapter.
To prepare the bed before the plate starts moving, particles are first created
and put into a container as in Chapter 2. The box shaped container’s cross
section is 40 cm $\times$ 30 cm and it is filled by particles to a depth of
approximately 30 cm. Once particles are created and placed, they are given
random initial velocities in three directions ($x$, $y$ and $z$). The absolute
velocity values for each direction are from 40 cm/s to 45 cm/s and the sign of
the velocity is randomly chosen. I wait 2.5 s in order to make sure the
granular bed is fully settled before starting the intrusion of the foot.
Figure 3.1 shows the mean speed values of 95700 particles and their standard
deviation from 0 to 9.9 s. The logarithmic mean and standard deviation can be
found in Figure 3.2. As the plots show, particle velocities quickly decay to
near 0, so 2.5 s is more than enough time to ensure the granular bed is fully
settled.

Figure 3.1: Average particle speed vs. time (upper) and particle speed
standard deviation vs. time (bottom). Figure 3.2: Logarithmic average particle
speed vs. time (upper) and logarithmic particle speed standard deviation vs.
time (bottom).

### 3.2 Resistive forces fitting functions

To get analytic expressions for $\alpha_{x,z}$ in Eq.2.1, I did hundreds of
DEM simulations with different attack angles $\beta$ and intrusion angles
$\gamma$ which vary from -90 degrees to 90 degrees. The constant intrusion
speed is chosen as 3 cm/s. For each $\beta-\gamma$ pair, I did three repeat
trials, and took the average of $\alpha_{x,z}$, as in Chapter 2, and then I
took absolute values of those $\alpha_{x,z}(\beta,\gamma)$ data. Figure 3.3
plots $|\alpha_{x,z}(\beta,\gamma)|$. In Figure 3.3, the upper two plots show
the $\alpha_{z}$ and $\alpha_{x}$ change with different $\beta$ and $\gamma$.
I took absolute values of the raw data because the stress results change
rapidly when the plate switches from intruding to extracting because the
stresses changes from positive to negative. For a small number of modes, the
Fourier transform is unable to capture this rapid change. However, since the
resistive forces are always opposite to the plate’s moving direction, I can
easily find the correct resistive forces direction given the plate’s
velocities. Thus, it’s reasonable to take absolute values of the raw data to
reduce Fourier fitting errors. By performing discrete Fourier transform on the
processed raw $\alpha_{z,x}$ data over $-\pi/2<\beta<\pi/2$ and
$-\pi<\gamma<\pi$, I can obtain a fitting function. To balance the accuracy
and conciseness of the fitting formulas, I set the Fourier transform order as
2. The fitting functions are:

$\alpha_{z}(\beta,\gamma)=\sum_{m=-2}^{2}\sum_{n=-2}^{2}[A_{m,n}\cos(2m\beta+n\gamma)+B_{m,n}\sin(2m\beta+n\gamma)]\mathrm{,}$
(3.1)
$\alpha_{x}(\beta,\gamma)=\sum_{m=-2}^{2}\sum_{n=-2}^{2}[C_{m,n}\cos(2m\beta+n\gamma)+D_{m,n}\sin(2m\beta+n\gamma)],$
(3.2)

where $A_{m,n}$, $B_{m,n}$, $C_{m,n}$ and $D_{m,n}$ are coefficients of the
fitting function and can be found in Table 3.2. The coefficients for 1st order
and 3rd order Fourier fitting can be found in this Google drive file
222https://drive.google.com/file/d/1T4itiIVEqeFl1zM6uPIz8_AqbE3kJiTy/view?usp=sharing.

$m$ | $n$ | $A_{m,n}$ | $B_{m,n}$ | $C_{m,n}$ | $D_{m,n}$
---|---|---|---|---|---
0 | 0 | 0.0587405 | 0 | 0.0510357 | 0
0 | 1 | -0.000449 | 0.0256221 | 0.0002929 | 0.0062092
0 | 2 | -0.007932 | -0.000897 | 0.0021059 | -0.0005
0 | -2 | -0.007932 | 0.0008966 | 0.0021059 | 0.0004998
0 | -1 | -0.000449 | -0.025622 | 0.0002929 | -0.006209
-1 | 0 | 0.0308869 | -0.0081 | -0.009617 | -0.007194
-1 | 1 | 0.0023856 | 0.0152946 | 0.0026648 | 0.0054052
-1 | 2 | -0.003031 | -0.005611 | -0.002222 | -0.003655
-1 | -2 | -0.00331 | -0.004783 | -0.002184 | -0.003597
-1 | -1 | -0.002843 | -0.016151 | -0.002984 | -0.005976
-2 | 0 | 0.0050729 | -0.00493 | -0.007099 | -0.002165
-2 | 1 | 0.0013374 | 0.0052393 | 0.0015156 | -0.005303
-2 | 2 | -1.17E-05 | -0.003304 | 0.0020161 | -0.000657
-2 | -2 | -0.000354 | -0.003454 | 0.0017459 | -0.001073
-2 | -1 | -0.001479 | -0.005812 | -0.001452 | 0.0051538
2 | 0 | 0.0050729 | 0.0049298 | -0.007099 | 0.0021647
2 | 1 | -0.001479 | 0.0058119 | -0.001452 | -0.005154
2 | 2 | -0.000354 | 0.0034541 | 0.0017459 | 0.0010725
2 | -2 | -1.17E-05 | 0.0033043 | 0.0020161 | 0.0006572
2 | -1 | 0.0013374 | -0.005239 | 0.0015156 | 0.0053032
1 | 0 | 0.0308869 | 0.0080999 | -0.009617 | 0.0071942
1 | 1 | -0.002843 | 0.0161505 | -0.002984 | 0.0059762
1 | 2 | -0.00331 | 0.0047832 | -0.002184 | 0.0035974
1 | -2 | -0.003031 | 0.0056105 | -0.002222 | 0.0036552
1 | -1 | 0.0023856 | -0.015295 | 0.0026648 | -0.005405
Table 3.2: Coefficients for 2nd order Fourier filtered force functions. Figure
3.3: Raw (top row) and Fourier filtered (bottom row) $\alpha_{z,x}$ colormaps
($\gamma$ and $\beta$ values are in radians). Colorbar units are in
N/$\textrm{cm}^{3}$.

Figure 3.4:  1st order $\alpha_{z,x}$ error maps ($\gamma$ and $\beta$ values
are in radians). Colorbar units are in N/$\textrm{cm}^{3}$. Figure 3.5: 2nd
order $\alpha_{z,x}$ error maps ($\gamma$ and $\beta$ values are in radians).
Colorbar units are in N/$\textrm{cm}^{3}$. Figure 3.6: 3rd order
$\alpha_{z,x}$ error maps ($\gamma$ and $\beta$ values are in radians).
Colorbar units are in N/$\textrm{cm}^{3}$. Figure 3.7: Mean error between
$\alpha_{x,z}$ from raw data and from Fourier representation vs. number of
modes.

To determine the accuracy of the fitting formulas, I generated error maps as
shown in Figure 3.4 to Figure 3.6. Errors for each $\beta-\gamma$ pair in the
error maps are computed as: $\alpha_{z}$ ($\alpha_{x}$) $-$ Fourier fitted
$\alpha_{z}$ (Fourier fitted $\alpha_{x}$). By comparing the error maps of
three Fourier fitting orders, one can see that a second order fit ensures the
fitting formula accuracy while limiting the number of fitting parameters. As
Figure 3.5 shows, errors in most regions are close to zero (less than 0.02)
and only errors of $\alpha_{z}$ near $\gamma=\pm\pi$ are near 0.1
$\mathrm{N}/\mathrm{cm}^{3}$. At the boundaries of the $\beta-\gamma$ plane,
the physical meaning is the plate switches between extraction and intrusion.
It is more difficult to intrude than withdraw, so the values near boundaries
see a rapid change along $\gamma$. However, low order Fourier fits are less
accurate if the raw data changes rapidly. This is why the fit is less accurate
near boundaries. Figure 3.7 shows that the mean and standard deviation of the
absolute errors for order 2 is less than 0.0211 $\mathrm{N/cm}^{3}$, so the
2nd-order Fourier series approximation can be used as a fast-to-evaluate model
of granular resistive stresses, at least for low intrusion velocities and away
from the effects of boundaries.

## CHAPTER 4Periodic Hopping of a Monopod

Many applications for mobile robots involve environments that are dangerous or
unsuitable for humans. Moreover, due to the inaccessible/unsafe nature of
these environments, level rigid ground is rare, so wheeled/treaded robots are
often unsuitable due to limited mobility. Legged robots would seem the ideal
choice, except that the bulk of research and development of legged robots has
focused on hard-ground applications, which are of little use in these
challenging environments. In many cases, the terrain is not rigid but, rather,
is some kind of deformable granular substrate (e.g., sand, soil, or snow).
While there has been recent progress on legged robotic locomotion on these
kinds of substrates (e.g., [3], [13], [5], [6]), a major remaining challenge
with hopping/walking/running on these substrates is that they do not return to
their pre-impact state. This represents a permanent energy loss as well as a
challenge to maintaining stable gaits. In addition to its relevance to legged
locomotion on deformable terrain, this unidirectional ground response
represents an interesting variation on the classic “bouncing ball problem” in
hybrid dynamics [14],[4]. While there are many tools for generating and
analyzing cyclic trajectories in hybrid dynamical systems (e.g., Poincare
analysis, Lyapunov analysis), I bring the tools of trajectory optimization [5]
and discrete-element-method (DEM) simulation to bear on the problem.
Specifically, I formulate the search for vertically-constrained hopping gaits
on deformable terrain as a trajectory optimization problem with boundary
constraints. I solve this problem numerically, obtaining a feedforward control
signal. Then, in DEM simulation, I apply this feedforward control signal along
with a stabilizing time-invariant feedback law.
In this chapter, I develop a control strategy for a vertically-constrained
monopod robot hopping on non-cohesive frictional granular media. The robot is
released from the apex of its trajectory and then impacts the bed with
velocities ranging from 0 to -4 m/s. When the foot contacts the granular bed,
the ground reaction forces are modeled by RFT. Then the equations of motion of
the robot has a closed-form solution. The trajectory optimization treats the
equations of motion as a part of the equality constraints and produces
feedforward control efforts for the robot. Finally, the center of mass (CoM)
of the robot can be controlled to return to the initial apex with zero speed
after both a single hop and five hops. Both single hop and five hops are
tested in DEM using both feedforward and feedforward + feedback control
strategies.

Figure 4.1: Physical models of the monopod robot and soft ground. The soft
ground is treated as a unidirectional spring with stiffness $k_{g}$, see Eq.
4.1. The body, located at height $q_{b}$ above the undisturbed ground, has
mass $m_{b}$, and the foot, located at height -$q_{f}$ below the undisturbed
ground, has mass $m_{f}$. Body and foot are connected with a linear motor
which exerts a possibly time dependent and equal and opposite force $u$ on
both masses. Heights $q_{b}$ and -$q_{f}$ are measured relative to the
undisturbed ground surface.

### 4.1 Soft ground model

Eq. 2.1 shows that the horizontal and vertical resistive forces increase
linearly with the intrusion depth. For a 1D hopping monopod robot, the
granular substrate can be modeled as a unidirectional spring. The ground
resistive force $f_{g}$ is expressed as

$f_{g}=\left\\{\begin{array}[]{ll}0&\textbf{flight}\\\
-k_{g}q_{f}&\textbf{yielding}\\\
\textbf{[}0,-k_{g}q_{f}\textbf{]}&\textbf{static}\end{array}\right.$ (4.1)

where $k_{g}$ is the ground stiffness. The area of the foot is 25
$\textrm{cm}^{2}$ and the stress gradient
$\alpha_{z}(\beta=0,\gamma=90^{\circ}$) is 0.25 $\textrm{N/cm}^{3}$. I chose
$k_{g}=25\textrm{ cm}^{2}\times 0.25\textrm{ N/cm}^{3}\times
100=625\mathrm{N/m}$ to match the experimental data in Chapter 3 I ignore the
withdraw force when lift off the robot. The foot thickness is chosen to ensure
no particles can accumulate on the top of the foot, which helps to reduce the
withdraw stresses.

### 4.2 Robot model

To focus on foot-ground interaction rather than whole-body control, I study a
simple monopedal robot hopping vertically as illustrated in Figure 4.1. The
robot consists of a body (position $q_{b}$ and mass $m_{b}$) and a flat-
bottomed foot (position $q_{f}$ and mass $m_{f}$) connected by a linear motor
(idealized as a source of force $u$) such that $u>0$ pushes the two masses
apart. The dimensions and masses of the robot foot and body can be found in
Table 4.1. The stroke limits of the robot system are

$0<q_{b}-q_{f}<0.5\ \mathrm{m.}$ (4.2)

Foot dimension | 5 cm x 5 cm x 10 cm
---|---
Foot mass | 0.25 kg
Body dimension | 5 cm x 5 cm x 5 cm
Body mass | 1.25 kg
Stroke limit | 0.5 m
Table 4.1: Physical properties of the monopod robot.

The hybrid dynamics of the hopping robot are divided into three phases:
flight, yielding stance, and static stance. Figure 4.2, reproduced from [8],
describes the conditions for transitioning between flight, yielding stance,
and static stance. Phase transitions depend on the foot kinematics and the
control effort. My robotic task is represented by periodicity constraints on
the CoM position and velocity, but during flight, I have no control over the
CoM trajectory. Thus, any attempt to influence the flight-phase CoM trajectory
must happen during stance. Therefore, I restrict my analysis to the static and
yielding stance phases. The equations of motion for the two stance phases are,

$\ddot{q_{b}}=\frac{u}{m_{b}}-g\mathrm{,and}$ (4.3)
$\ddot{q_{f}}=\frac{f_{g}}{m_{f}}-\frac{u}{m_{f}}-g\mathrm{,}$ (4.4)

where $g$ is the gravitational acceleration ($9.81\
\mathrm{m}/\mathrm{s}^{2}$). For the following section on robot control,
initial conditions are,

$\begin{array}[]{c}q_{b}(0)=0.25\textrm{ m, }q_{f}(0)=0\textrm{ and
}\dot{q}_{b}=\dot{q}_{f}=V_{0}\textrm{,}\end{array}$ (4.5)

where $V_{0}$ is the impact velocity.

Figure 4.2: Transitions between flight phase, yielding stance and static
stance. Reproduced from [8] with permission.

### 4.3 Robot control

My first task is to find period one hopping gaits for the vertically-
constrained monopod robot. Gait generation is the formulation and selection of
a sequence of coordinated leg and body motions that propel a legged robot
along a desired path [15]. Specifically, my robot starts impact at a certain
initial position and moves vertically downward with an initial impact
velocity. The desired motor control periodically returns the robot to the
initial impact position with the same CoM velocity magnitude but in the upward
direction on each successive hop.
There are various methods of determining the control signal $u(t)$ that
results in periodic hopping trajectories (e.g., Poincare analysis [16],
Lyapunov analysis), but I use trajectory optimization [17] to find controls
and trajectories that satisfy periodicity constraints.

####  4.3.1Feedforwrd control

To determine the open loop control signal for the robot, I divide the problem
into two parts, stance phase and flight phase. The stance-phase optimal
control problem can be transformed into a constrained nonlinear program. For
the flight phase, no ground forces are applied to the foot, the velocity of
body and foot is expected to be equal, and the distance between the body and
foot should return to the initial value, so analytic methods are used to
determine the corresponding control signal.
Hubicki et al. [5] describe a fast optimal motion planning algorithm for a 1D
hopping robot. By formulating an optimal control problem as a constrained
nonlinear program (NLP), I utilize a well-developed NLP solver (i.e., MATLAB
fmincon) to find out the control signal. Control signals $u(t)$ are chosen as
fourth order polynomials because the system has 4 constraints, which determine
4 coefficients of the polynomials. The remaining coefficient is determined by
optimizing the cost function. The expression is as follows,

$u(t)=a+bt+ct^{2}+dt^{3}+et^{4}\mathrm{,}$ (4.6)

where $\mathbf{x}=[a,b,c,d,e]$ are the design variables of the NLP. Note that
I fix the time at 0.3 s for the stance phase and optimize the coefficients of
the control signals $u(t)$ to satisfy the constraints. I choose 0.3 s because
it can ensure my MATLAB code to find a solution for the optimal control. For a
longer time like 1s can also make the optimization problem converge but it
needs more time to solve. In addition, too short a stance time will not be
able to develop enough momentum for the robot to jump back to apex given
stroke constraints and lack of rate dependence in GRF model.
I choose control effort squared as the cost function:

$J(\mathbf{x})=\int_{0}^{T_{s}}u^{2}(t)dt\textrm{,}$ (4.7)

Reference [8] chooses intrusion depth of the foot as the cost function. I
choose the control effort squared because it makes the cost function convex
which is good for the convergence of the optimization problem. The resulting
formation for the stance-phase nonlinear program is as follows:

$\begin{array}[]{lll}\mathbf{x}^{*}&=\displaystyle\operatorname*{arg\,min}_{\mathbf{x}}J(\mathbf{x})&\textrm{s.t.}\\\
f(\mathbf{x})&=0\textrm{,}\\\ g(\mathbf{x})&\leq 0\mathrm{,}\end{array}$ (4.8)

where $f(\mathbf{x})$ are the equality constraints obtained by the robot
dynamics (see Eq. 4.3 and Eq. 4.4 ) and robot boundary conditions,
$g(\mathbf{x})$ are the inequality constraints obtained by the linear
controller control limits and robot stroke limits (see 4.2), and $T_{s}$ is
the desired stance-phase period.
The stance-phase terminal conditions at $T_{s}$ are

$\begin{array}[]{cc}q_{b}(T_{s})=C\textrm{, }q_{f}(T_{s})=D\textrm{,
}\dot{q}_{b}(T_{s})=E\textrm{ and }\dot{q}_{f}(T_{s})=0\mathrm{,}\end{array}$
(4.9)

where $[C$, $D$, $E]$ (terminal body position, terminal foot position, and
terminal body velocity respectively) are constants determined by ground
stiffness and impact velocity. One thing interesting is my robot never leaves
ground in the stance phase because I add constraints in the MATLAB solver to
do this.
The flight phase control first applies a constant forces to the foot and body,
and then turns off the linear motor during the remainder of free fall. I
choose this strategy for the flight phase because the solution for control
efforts is easily obtained. Once given the depth of the foot at the end of the
stance phase and the lift-off acceleration of the foot, I can solve for the
control efforts for lift-off and the terminal states of the body in the stance
phase.
By setting the stance-phase duration to 0.3 s and boundary conditions as given
in Table 4.2, I obtain the control efforts and the resulting body and foot
trajectories shown in Figure 4.3. As those results are for an example impact
velocity of -0.2 m/s, the foot and the body move downwards from the top of the
granular media. The control signals make the foot and the body return to the
initial height while the terminal speed values are the same as the initial
speed. We can also observe the double intrusion of the foot, which is an
interesting feature of locomoting on soft ground.

Figure 4.3: Planned foot and body states vs. time (upper and middle plots) and
planned control efforts vs. time (bottom plot).

$q_{f}(0)$ | 0 m | $\dot{q}_{f}(0)$ | -0.2 m/s
---|---|---|---
$q_{b}(0)$ | 0.25 m | $\dot{q}_{b}(0)$ | -0.2 m/s
$q_{f}(T)$ | 0 m | $\dot{q}_{f}(T)$ | 0.2 m/s
$q_{b}(T)$ | 0.25 m | $\dot{q}_{b}(T)$ | 0.2 m/s
Table 4.2: Boundary conditions

####  4.3.2Feedforward plus feedback control

After testing the open loop control forces in Chrono simulation, I found that
the trajectories of the robot deviate from the expected ones because the
ground stiffness approximation is not sufficiently accurate. Therefore, I
introduce a feedback controller to stabilize control in stance phase. So, the
control strategy is feedforward plus feedback control now.
As Figure 4.4 shows, the control system inputs are the desired center of
mass(CoM) position and velocity of the robot. The desired CoM
states($q_{CoM,d}\textrm{ and }\dot{q}_{CoM,d}$) are generated from MATLAB
when I apply the open loop control forces to the robot in numerical
simulation. The difference between the desired robot states and the actual
robot states ($q_{CoM}\textrm{ and }\dot{q}_{CoM}$) are the CoM
errors($q_{CoM,e}\textrm{ and }\dot{q}_{CoM,e}$). By passing those errors into
a PD controller, I obtain the feedback control forces $U_{FB}$ as illustrated
in Fig.4.4. In the section above, I use trajectory optimization to get the
feedforward control forces($U_{FF}$). Combining the feedback and feedforward
control forces and adding the new forces($U$) into the equations of motion of
the robot system, I get the real-time states of the robot.
I implement this feedforward plus feedback controller on my robot model in
Chrono. Specifically, I create a robot model like Figure 4.1 with physical
properties as given in Table 4.1. I set the impact velocity of both the body
and the foot to -0.2 m/s. When the robot transitions from stance phase to
flight phase, the CoM position is expected to return to 0, and the CoM
velocity is expected to be 0.2 m/s. Since no external force is applied to the
robot system during the flight phase, the robot can return to the initial
position with the same velocity because of momentum conservation. The two
conditions thus ensure the robot hops periodically. For the open loop control
case, I applied the feedforward control forces generated from MATLAB. For the
feedforward plus feedback case, I add a PD controller on the robot to track
the CoM of the robot. As Figure 4.5 shows, my controller significantly reduce
both the CoM position and velocity errors.
To quantitatively analyze the error of the two control strategies (i.e.
feedforward and feedforward+feedback), I test their performance for different
impact velocities(from -4 m/s to -0). The errors are computed as the terminal
CoM position/velocity of the two control methods in Chrono minus the desired
terminal CoM position/velocity from MATLAB. For the feedforward plus feedback
control, I run Chrono simulation for 3 trials to analyze its uncertainty. As
Figure 4.6 shows, my strategy works well in reducing the error for different
impact velocities compared to the open loop control alone.
In addition, I investigate how different foot sizes influence the performance
of my controllers. As the foot size decreases toward the particle size, force
fluctuations increase and present greater challenges to the controllers. To
investigate the efficacy of the controller, I test three foot sizes (1 cm x 1
cm, etc.) in addition to the original 5 cm x 5 cm foot. To ensure the foot
reach the same depth for different foot sizes, I scale the mass of the foot
and the body based on the nondimensionalization method in [8]. The ratio of
the total mass of the robot to the foot square area is fixed. For each foot
size, I test six different impact velocities(0, -0.2m/s, -0.5m/s, -1m/s, -2m/s
and -4m/s) and conduct three trials for each impact velocity. Both the
feedforward + feedback controller and the open loop controller are tested in
these experiments. The CoM position error percentage is computed as Eqs. 4.10,

$\begin{array}[]{ll}\eta&=\frac{q_{CoM,i}(T)-q_{CoM,d}(T)}{q_{CoM,d}(T)}\times
100\%,\end{array}$ (4.10)

where $i$ can be feedforward + feedback or feedforward alone, $d$ is the
desired CoM position generated from MATLAB, and $T$ is the duration of one
hop. For each foot size, I compute the mean values and standard deviations of
all the error percentages of all impact velocities, and draw errorbars as
shown in Figure 4.7. Again, I observe that the feedforward plus feedback
controller reduces the error to less than 2$\%$ in all cases. This indicates
my controller is robust for different foot sizes, even in the presence of
large variations in the ground reaction force when the foot is small.

Figure 4.4: Block diagram of feedforward plus feedback control for stance
phase. $q_{CoM}$ and $\dot{q}_{CoM}$ are the center of mass position and
velocity, respectively. $q_{CoM}$, $q_{CoM,d}$, and $q_{CoM,e}$ are actual
center of mass position, desired position and position error, respectively.
$U_{FB}$ and $U_{FF}$ are the feedback and feedforward control forces,
respectively. Figure 4.5: CoM trajectories from optimal(desired) and for two
control strategies. Upper plot: CoM position vs. time. Bottom plot: CoM
velocity vs. time. Blue, red and yellow curves are feedforward plus feedback
results, expected results from MATLAB and feedforward control results,
respectively. Figure 4.6: Position and velocity errors of two different
control strategies for multiple impact velocities. Asterisks are errors of
open loop control while error bars are errors of feedforward plus feedback
control. Figure 4.7: Experimental errors in CoM position from Eq.4.10 vs. foot
size.

### 4.4 Multiple hops

The ultimate goal of my work is to design gaits for monopod hopping on
deformable ground. For the single hop in the sections above, the robot jumps
on undisturbed ground. To challenge my feedforward + feedback controller, I
make the monopod to have five hops on disturbed ground. As Section 4.3.2
shows, I successfully implemented my feedforward plus feedback controller on
the robot for single hop. Figure 4.3 shows expected trajectories of a single
hop of the monopod. The robot starts hopping with impact velocities of -0.2
m/s and returns to the initial position with the same velocity. My control
strategy for mutiple hops is to track the CoM trajectories of the single hop
of the robot in each hop. In each hop of a multi-hop sequence, I set the
period equal to the duration of expected trajectories from MATLAB and use my
feedforward plus feedback controller to stabilize my control in stance phase.
For the flight phase, I add a PD controller to make the velocity of the body
and foot equal, and make the body-foot distance equal to the initial value.
Figure 4.8 compares the CoM trajectories of two control strategies. The robot
starts moving with an impact velocity of -0.2 m/s. The feedforward plus
feedback controller has far greater performance than the open loop controller.
After five hops, the CoM trajectories of the robot tracks the expected
trajectories quite well.
To compare the errors associated with the two control strategies, I compute
the errors as CoM position/velocity at the end of five hops minus the expected
terminal CoM position/velocity. To ensure my controller is robust for
different impact velocities, I test impact velocities ranging from -4 m/s to
0. For each impact velocity, I run Chrono simulation for three trials. Figure
4.9 shows errorbars of both CoM position and CoM velocity after five hops for
each impact velocity. The errors of the feedforward plus feedback controller
are quite close to zero for all the impact velocities I tested. This again
shows the robustness of my controller. In addition to the plots in this
chapter, I also made corresponding videos of the robot hopping in Chrono. The
first
video111https://drive.google.com/file/d/1LkUWFWx_U1X0XtGNpdf6Dcj_IhaWzxXH/view?usp=sharing
shows five hops on soft ground using a feedforward plus feedback controller.
The robot in the second
video222https://drive.google.com/file/d/1E8iMqUuT_bcmHaP9IYCvcwEjnf4R89iO/view?usp=sharing
is controlled by an open loop controller. Figure 4.10 is the screenshot of one
demo video. Comparing the two videos shows that the feedback controller works
well making the CoM to reach the expected height.

Figure 4.8: Center of mass trajectories for five hops. Upper plot: CoM
position vs. time. Bottom plot: CoM velocity vs. time. Blue, red and yellow
curves are feedforward plus feedback results from Chrono, expected results
from MATLAB, and feedforward control results from Chrono, respectively. Figure
4.9: Errorbar of robot CoM trajectories after five hops for different impact
velocities under open loop (orange) and feedforward plus feedback (blue)
control. Figure 4.10: Screenshot of demo video. Upper yellow box is the body
while bottom yellow box is the foot. The blue line represents the linear motor
which connects the body and the foot. Red sphere is the real time position of
the CoM of the robot. The blue sphere is the desired apex for the CoM.

## CHAPTER 5Conclusion

In this work, I set up a DEM based simulation environment called Chrono. To
validate its reliability and accuracy, I ran many simulations and compared my
results with results in [3]. The simulation results have a good match to
experimental results. After that, I utilized Chrono to generate a ground
stiffness approximation formula based the ground resistive force model in [3].
I also formulated a soft-landing problem and develop solutions to design open
loop periodic gait for our monopod hopper with the trajectory optimization and
the approximation formula. Then I develop a new control solution by applying
feedforward plus feedback control on the monopod robot. The new control
strategy is tested on single hops and multiple hops for various impact
velocities in Chrono simulation. The Chrono simulation results indicate our
controller is robust and can reduce CoM trajectory errors.

This work proves controlling a monopod to desired height periodically when
hopping on yielding terrain is possible. However, our solution for the control
force is not optimal, at least in terms of the energy lost to ground. My
feature plan is to seek the optimal control that ensures the periodically
hopping while minimizes the intrusion depth. Besides, all the experimental
results I got are based on DEM simulation. It’s interesting to carry out
experiment validation on a real robot.

## Bibliography

  * [1] J. W. Grizzle et al.  “MABEL, a new robotic bipedal walker and runner”  In _2009 American Control Conference_ , 2009, pp. 2030–2036  DOI: 10.1109/ACC.2009.5160550
  * [2] Ian R. Manchester, Uwe Mettin, Fumiya Iida and Russ Tedrake  “Stable Dynamic Walking over Rough Terrain”  In _Robotics Research_ Berlin, Heidelberg: Springer Berlin Heidelberg, 2011, pp. 123–138
  * [3] Chen Li, Tingnan Zhang and Daniel I. Goldman  “A Terradynamics of Legged Locomotion on Granular Media”  In _Science_ 339.6126 American Association for the Advancement of Science, 2013, pp. 1408–1412
  * [4] Jeffrey Aguilar, Alex Lesov, Kurt Wiesenfeld and Daniel I. Goldman  “Lift-Off Dynamics in a Simple Jumping Robot”  In _Phys. Rev. Lett._ 109 American Physical Society, 2012, pp. 174301  DOI: 10.1103/PhysRevLett.109.174301
  * [5] C. M. Hubicki, J. J. Aguilar, D. I. Goldman and A. D. Ames  “Tractable terrain-aware motion planning on granular media: An impulsive jumping study”  In _2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2016, pp. 3887–3892  DOI: 10.1109/IROS.2016.7759572
  * [6] X. Xiong, A. D. Ames and D. I. Goldman  “A stability region criterion for flat-footed bipedal walking on deformable granular terrain”  In _2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2017, pp. 4552–4559  DOI: 10.1109/IROS.2017.8206323
  * [7] Sonia F Roberts and Daniel E Koditschek  “Reactive Velocity Control Reduces Energetic Cost of Jumping with a Virtual Leg Spring on Simulated Granular Media”  In _Proceedings of the 2018 IEEE International Conference on Robotics and Biomimetics (ROBIO)_ , 2018  IEEE
  * [8] D. J. Lynch, K. M. Lynch and P. B. Umbanhowar  “The Soft-Landing Problem: Minimizing Energy Loss by a Legged Robot Impacting Yielding Terrain”  In _IEEE Robotics and Automation Letters_ 5.2, 2020, pp. 3658–3665  DOI: 10.1109/LRA.2020.2977260
  * [9] Matthias Sperl  “Experiments on Corn Pressure in Silo Cells – Translation and Comment of Janssen’s Paper from 1895”  In _Granular Matter_ 8, 2005  DOI: 10.1007/s10035-005-0224-z
  * [10] Shinichiro Miyai, Murino Kobayakawa, Takuya Tsuji and Toshitsugu Tanaka  “Influence of particle size on vertical plate penetration into dense cohesionless granular materials (large-scale DEM simulation using real particle size)”  In _Granular Matter_ 21, 2019  DOI: 10.1007/s10035-019-0961-z
  * [11] Paul Umbanhowar and Daniel I. Goldman  “Granular impact and the critical packing state”  In _Phys. Rev. E_ 82 American Physical Society, 2010, pp. 010301  DOI: 10.1103/PhysRevE.82.010301
  * [12] Feifei Qian et al.  “Principles of appendage design in robots and animals determining terradynamic performance on flowable ground”  In _Bioinspiration & Biomimetics_ 10.5 IOP Publishing, 2015, pp. 056014  DOI: 10.1088/1748-3190/10/5/056014
  * [13] J. Aguilar and D. Goldman  “Robophysical study of jumping dynamics on granular media”  In _Nature Physics_ 12, 2016, pp. 278–283
  * [14] P.J. Holmes  “The dynamics of repeated impacts with a sinusoidally vibrating table”  In _Journal of Sound and Vibration_ 84.2, 1982, pp. 173–189  DOI: https://doi.org/10.1016/S0022-460X(82)80002-3
  * [15] D. Wettergreen and C. Thorpe  “Gait Generation For Legged Robots”  In _Proceedings of the IEEE/RSJ International Conference on Intelligent Robots and Systems_ 2, 1992, pp. 1413–1420  DOI: 10.1109/IROS.1992.594568
  * [16] Eric R Westervelt et al.  “Feedback control of dynamic bipedal robot locomotion”  CRC press, 2018
  * [17] Matthew Kelly  “An Introduction to Trajectory Optimization: How to Do Your Own Direct Collocation”  In _SIAM Review_ 59.4, 2017, pp. 849–904


# ILMART: Interpretable Ranking with Constrained LambdaMART

Claudio Lucchese<EMAIL_ADDRESS>0000-0002-2545-0425 Ca’ Foscari
University of VeniceItaly ,  Franco Maria Nardini
<EMAIL_ADDRESS>0000-0003-3183-334X ISTI-CNRItaly ,  Salvatore
Orlando<EMAIL_ADDRESS>0000-0002-4155-9797 Ca’ Foscari University of
VeniceItaly ,  Raffaele Perego  0000-0001-7189-4724
<EMAIL_ADDRESS>ISTI-CNRItaly and  Alberto Veneri
<EMAIL_ADDRESS>0000-0003-2094-3375 Ca’ Foscari University of
VeniceISTI-CNRItaly

###### Abstract.

Interpretable Learning to Rank (LtR) is an emerging field within the research
area of explainable AI, aiming at developing intelligible and accurate
predictive models. While most of the previous research efforts focus on
creating post-hoc explanations, in this paper we investigate how to train
effective and intrinsically-interpretable ranking models. Developing these
models is particularly challenging and it also requires finding a trade-off
between ranking quality and model complexity. State-of-the-art rankers, made
of either large ensembles of trees or several neural layers, exploit in fact
an unlimited number of feature interactions making them black boxes. Previous
approaches on intrinsically-interpretable ranking models address this issue by
avoiding interactions between features thus paying a significant performance
drop with respect to full-complexity models. Conversely, ILMART, our novel and
interpretable LtR solution based on LambdaMART, is able to train effective and
intelligible models by exploiting a limited and controlled number of pairwise
feature interactions. Exhaustive and reproducible experiments conducted on
three publicly-available LtR datasets show that ILMART outperforms the current
state-of-the-art solution for interpretable ranking of a large margin with a
gain of nDCG of up to 8%.

Interpretable Ranking, Interpretable Boosting, LambdaMART

††booktitle: Woodstock ’18: ACM Symposium on Neural Gaze Detection, June
03–05, 2018, Woodstock, NY††price: 15.00††isbn: 978-1-4503-XXXX-X/18/06††ccs:
Information systems Learning to rank††ccs: Computing methodologies Boosting

##  1\. Introduction

Learning to Rank (LtR) is a vast research area focusing on learning effective
models for solving ranking tasks. Although being effective, state-of-the-art
LtR techniques do not deal with the _interpretability_ of the learned model,
i.e., the possibility given to humans to easily understand the reasoning
behind predictions made by the model. The need for learning trustworthy models
in LtR is becoming a topic of increasing interest, which tries to address
important concerns about privacy, transparency, and fairness of ranking
processes behind many everyday human activities.

For what regards the fairness of ranking algorithms, several solutions have
been proposed to create fairer models, where the term fair takes different
meanings according to the metric employed (Castillo, 2019). However, focusing
only on devising a fairer model does not imply to be compliant with new
regulations that are emerging. An example of that is _The Right of
Explanation_ , which allows each European citizen to claim a justification for
every algorithm decision made with their data that affects their life (Politou
et al., 2018). Therefore, for highly-sensitive scenarios new LtR approaches
have to be investigated to train ranking models that do not jeopardize
interpretability in the name of effectiveness.

In literature, two main lines of work dealing with the interpretation of
ranking models are investigated. The first concerns the creation of post-hoc
explanations for black-box models, while the second focuses on the design of
intrinsically-interpretable models. In the first line, the goal is the
definition of techniques that allow, given a ranking model, to generate post-
hoc explanations of the decisions made. Here, some of the explanation
techniques used in regression or classification tasks have been adapted to the
ranking scenario. For example, two recent contributions (Singh and Anand,
2019; Verma and Ganguly, 2019) present two different approaches to adapt LIME
(Ribeiro et al., 2016) to the ranking task. The work by Singhs and Anand
focuses only on pointwise LtR techniques and treats the ranking problem as a
classification problem, where the probability of a document to be relevant is
computed with respect to its score, and a simple linear SVM is used as a local
explanation model (Singh and Anand, 2019). The work by Verma and Ganguly
evaluates the consistency and correctness of the explanations generated by
LIME by varying the sampling strategies used to create the local neighborhood
of a document (Verma and Ganguly, 2019). In the same line, Fernando _et al._
propose an approach to reuse SHAP (Lundberg and Lee, 2017) for the ranking
task (Fernando et al., 2019). Even though post-hoc explanations can be seen as
a viable way to inspect and audit accurate ranking models, they can produce
explanations that are not faithful to internal computations of the model, and
thus they can be misleading (Rudin, 2019). For this reason, Zhuang _et al._
pioneer the problem of developing intrinsically-interpretable ranking models
(Zhuang et al., 2021). In detail, they propose Neural RankGAM, an adaptation
of GAM (Hastie and Tibshirani, 1986) to the ranking task. Although Neural
RankGAM nicely exploits the interpretability of GAMs, the authors show that
the ranking performance achieved is significantly worse than those provided by
state-of-the-art black-box models. However, it is worth noticing that in other
contexts, such as in regression and classification tasks, other intrinsically
interpretable models appear to be as accurate as other black-box models. A
notable example of an interpretable and effective solution is Explainable
Boosting Machine (EBM), a C++ implementation of the $GA^{2}M$ algorithm (Lou
et al., 2013), where each component of the Generalized Additive Models (GAM)
are learned through a mixed bagging-boosting procedure (Lou et al., 2012).

Figure 1. ILMART learning process. Each color corresponds to a different
feature used in the splits, and nodes colored in black are leaves. The trees
for the main and interaction effects are represented with the same depth for
illustration purposes only.

In this paper, we contribute to the development of intrinsically-interpretable
ranking models with ILMART, a novel LtR solution based on LambdaMART. ILMART
creates effective and intelligible models by encompassing a limited and
controlled number of pairwise feature interactions. We evaluate ILMART on
three publicly-available LtR datasets. The results of our reproducible
experiments show that ILMART outperforms the current state-of-the-art solution
for interpretable ranking of a large margin, with a gain in terms of nDCG of
up to 8% with respect to Neural RankGAM.

##  2\. Interpretable LambdaMART

ILMART produces ensemble-based ranking models having the structure of a GAM,
which is considered an interpretable model (Lou et al., 2012). Before
introducing ILMART, we review the main aspects of GAMs to highlight
similarities and differences with respect to the newly proposed technique.

Generalized Additive Models. Let $\bm{x}\in\mathbb{R}^{d}$ be an instance in a
$d$-dimensional vector space, with $x_{j}$ being the value of its $j$-th
feature. An instance $\bm{x}$ is associated with a target label $y$ generated
by some unknown target function $g(\cdot)$, i.e., $g(\bm{x})=y$. In its basic
form, a GAM models the target variable $y$ as a combination of $p$ main
effects and $K$ non-zero interaction effects, as defined in (Hastie and
Tibshirani, 1986):

$l(\mu(\bm{x}))=\alpha+\sum_{j=1}^{p}s_{j}(x_{j})+\sum_{j=1,k=1}^{j=p,k=p}s_{jk}(x_{j},x_{k})$

Where $\mu(\bm{x})=E[y|\bm{x}]$, with $y$ following a distribution of the
exponential family, and $l(\cdot)$ being the so-called _link function_ that is
able to describe the relationship between the $p$ main effects and the $K$
interactions effects with the expected value of $y$. In addition, the number
of non-zero interaction effects is equal to $K$, i.e.
$|\\{s_{jk}(\cdot)|s_{jk}(\cdot)\neq 0\\}|=K$. Usually, the $p$ main effects
and the and the $K$ interactions are modeled through simple third order
splines. From the explainability point of view, the analyst is promptly
provided with a plot of each basis function $s_{j}$ that clearly describes the
relation between each feature and the target variable $y$, while $s_{jk}$
values represent the contribution of the interactions between feature pairs.
For classification and regression tasks, (shallow) decision trees instead of
splines were investigated in (Lou et al., 2013).

Interpretable LambdaMART. We now introduce ILMART, our novel method based on
LambdaMART to learn ranking models that are both accurate and interpretable.
ILMART presents similarities with GAM, and builds the prediction $\hat{y}$ of
a unknown target variable $y$ as follows:

(1)
$\hat{y}=\overbrace{\underbrace{\sum_{j\in{\mathcal{J}}}\tau_{j}(x_{j})}_{R\
\textrm{trees}}}^{|{\mathcal{J}}|=p\ \textrm{main
effects}}+\overbrace{\underbrace{\sum_{(i,j)\in{\mathcal{K}}}\tau_{ij}(x_{i},x_{j})}_{T\
\textrm{trees}}}^{|{\mathcal{K}}|=K\ \textrm{interaction effects}}$

where $\tau_{j}(x_{j})$ and $\tau_{ij}(x_{i},x_{j})$ are ensembles of trees.
Each tree $t\in\tau_{j}(x_{j})$ models a _main effect_ in which the only
feature allowed for each split is $x_{j}$.

On the other hand, each tree $t\in\tau_{ij}(x_{i},x_{j})$ models the
_interaction effects_ , where the only features allowed in a split of tree $t$
are either $x_{i}$ or $x_{j}$. Similarly to the notation used to describe a
GAM, $p$ denotes the number of distinct main effects $x_{i}$, whereas $K$
denotes the number of distinct interaction effects $(x_{i},x_{j})$.

To learn the ensemble of trees modeled by Eq. (1), we adopt a modified version
of the LambdaMART boosting algorithm. Unlike GAM, we do not make any
assumption on the distribution of $y$ to guide the learning process. The
learning strategy to generate the main effects $t\in\tau_{j}(x_{j})$ and the
interaction effects $t\in\tau_{ij}(x_{i},x_{j})$ in Equation 1 is made of
three steps depicted in Fig. 1: _i)_ Main Effects Learning, which learns a set
of trees, each one working on a main effect (single feature); _ii)_
Interaction Effects Selection, which selects the top-$K$ most important
interactions effects (feature pairs); _iii)_ Interaction Effects Learning,
which learns a set of trees, each one exploiting one of the interaction
effects selected by the previous step.

_Main Effects Learning_. In the initial step, illustrated on top of Fig. 1, we
learn an ensemble of $R$ trees modeling the main effects. To this end, we
constrain the LambdaMART boosting procedure to use a single feature per tree
in the first $R$ boosting rounds. In other words, at a given boosting round of
LambdaMART, when feature $x_{j}$ is chosen for the root’s split of a tree $t$
, we force the algorithm to use $x_{j}$ for all the further splits in $t$,
until a stopping criterion for tree growth is reached. LambdaMART stops its
boosting rounds when no improvements are observed on the validation set.
Eventually, we obtain an ensemble of $R$ trees only using for their splits a
set ${\mathcal{J}}$ of $p$ features, $p\leq d$ and $p\leq R$, modeling the
main effects of Equation 1. Thus, $p$ is not an input parameter of ILMART, but
it is the number of distinct features used by the ensemble of $R$ trees
modeling the main effects only. We name this initial model ILMART. It is worth
noting that, as a side effect of this first step, we achieve the additional
effect of potentially reducing the cardinality of the features, from $d$ to
$p$.

_Interaction Effects Selection_. In the second step, we still exploit the
LambdaMART boosting procedure to select the top-$K$ interaction effects, as
illustrated in the middle of Fig. 1. Specifically, we continue the learning of
the ILMART model obtained so far, by enabling the interactions between all
possible pairs of $p$ features identified by the previous step, according to
the interaction heredity principle (Cox, 1984). This is accomplished by
constraining LambdaMART to fit trees with only $3$ leaves (and $2$ splits),
where the $2$ splits use a pair of distinct features. The boosting round is
stopped when the size of the set ${\mathcal{K}}$, composed of the distinct
feature pairs used by the ensemble of $S$ trees, is exactly $K$, where
$K\leq\binom{p}{2}$. The $S$ trees are generated solely for the purpose of
identifying the interaction effects and they are discarded afterwards.

Note that for selecting the main $K$ interaction effects, we rely on the
capability of LambdaMART in identifying, at each boosting round, the features
and the splits with the highest contributions in minimizing the loss function.

Thus, the interaction effects selected by the boosting procedure are likely be
generated in order of importance.

_Interaction effects learning_. Finally, the last step sketched on the bottom
of Fig. 1 consists in adding to the initial ILMART model $T$ new trees learned
by constraining LambdaMART to use only the top-$K$ feature pairs identified in
the previous step. The number $T$ of new learned trees is not an input of
algorithm, but it is chosen on the basis of the validation set. The final
model, obtained by adding to ILMART the $T$ trees working on the $K$
interaction effects is named ILMARTi.

##  3\. Experimental Evaluation

We experimentally assess the performance of ILMART and ILMARTi on three public
datasets for LtR, namely Istella-S, Web30K, and Yahoo. The Istella-S dataset
(Lucchese et al., 2016) includes $33$,$018$ queries with an average of $103$
documents per query. Each document-query pair is represented by $220$
features. The Web30K dataset (Qin and Liu, 2013) is composed of more than
$30$,$000$ queries (in $5$ folds), with an average of $120$ documents per
query and $136$ features per query-document pair. The Yahoo dataset (Chapelle
and Chang, 2011) is composed of two sets, including document-query pairs with
$699$ features. The feature vectors of the three datasets are labeled with
relevance judgments ranging from $0$ (irrelevant) to $4$ (perfectly relevant).
In our experiments, we use train/validation and test splits from fold 1 of the
Web30K dataset and data from “set 1” of the Yahoo dataset.

Competitors. We compare the performance of ILMART and ILMARTi with the
following competitors:

  * •

Neural RankGAM (NRGAM), a neural-based approach to learn GAMs for ranking. The
approach exploits standalone neural networks to instantiate sub-models for
each individual feature. It is the current state-of-the-art technique for
learning interpretable models for ranking (Zhuang et al., 2021).

  * •

EBM, an efficient implementation of $GA^{2}M$ using a mixed bagging-boosting
procedure (Lou et al., 2013). In this case the learned model is used as a
pointwise ranker, since the implementations available are only made to solve
regression and classification tasks. We denote with EBMi the EBM model using
interaction effects.

Metrics. We measure the performance in terms of normalized Discounted
Cumulative Gain (nDCG) at three different cutoffs, i.e., $\\{1,5,10\\}$. We
compute the nDCG metric by employing exponential weighing of the relevance
(Burges et al., 2005). By default, queries with missing relevant documents
have been assigned a nDCG $=1.0$. The statistical significance is computed
with a two-sided Fisher’s randomization test (Smucker et al., 2007) with
significance level $p<0.05$. The test is computed using the RankEval library
(Lucchese et al., 2017).

Implementation and training settings. We implement ILMART as an extension of
the LightGBM library111https://github.com/microsoft/LightGBM (Ke et al.,
2017). The extension works by adding the possibility to constraint the
boosting procedure to use only a limited number of features per tree. The
source code used in our experiments along with the models trained on the three
public datasets is publicly available
online222https://github.com/veneres/ilmart.

The training of ILMART and ILMARTi optimizes nDCG@10. We perform hyper-
parameter tuning by varying the number of leaves in the trees in
$\\{32,64,128\\}$ and the learning rate in $\\{0.001,0.01,0.1\\}$. Early
stopping, i.e., stopping boosting if no improvement on the validation set is
observed for $100$ consecutive rounds, is used for training in step one (main
effects learning) and three (interaction effects learning). The optimal ILMART
models learned with the settings above are made of $914$, $895$ and $1$,$969$
trees for the Web30K, Yahoo, and the Istella-S datasets, respectively. On the
same datasets, the three optimal ILMARTi models are made of $1$,$365$,
$1$,$215$ and $3$,$024$ trees, respectively. Neural RankGAM is also trained by
optimizing nDCG@10. We use the public implementation of NRGAM available in
Keras333https://github.com/tensorflow/ranking/issues/202 and made available by
the authors of the original paper (Zhuang et al., 2021). We train NRGAM by
using the same settings reported by the authors. In addition, we employ a
batch size of $128$, and we train the networks with $3$,$000$ epochs for
Web30K, $700$ for Yahoo and $1$,$000$ for Istella-S, with early stopping set
to $100$ epochs. We also learn EBM by optimizing nDCG@10. The performance of
EBM is fine-tuned by varying the number of outer-bags in {5, 10, 15} for the
two versions with and without interactions. In our experiments, we employ a
publicly-available implementation of EBM
444https://interpret.ml/docs/ebm.html. For what regards the two strategies
encompassing interaction effects, i.e., ILMARTi and EBMi, we limit the number
of maximum interactions allowed to $50$.

Table 1. Performance comparison in terms of nDCG. Statistically significant improvements w.r.t. NRGAM are marked with an asterisk (*). Best results are reported in bold. Dataset | Method | nDCG | $p$ | $K$
---|---|---|---|---
@1 | @5 | @10
Web30K | NRGAM | 44.80 | 43.72 | 45.73 | 136 | 0
EBM | 41.70 | 43.14 | 45.66 | 136 | 0
EBMi | 46.51* | 45.96* | 48.01* | 136 | 50
ILMART | 45.17 | 45.00* | 47.05* | 79 | 0
ILMARTi |  48.94* |  47.76* |  49.55* | 79 | 46
Yahoo | NRGAM | 67.92 | 70.13 | 75.07 | 699 | 0
EBM | 70.33* | 72.77* | 77.29* | 699 | 0
EBMi | 70.22* | 72.85* | 77.36* | 699 | 50
ILMART | 70.57* | 72.64* | 77.21* | 153 | 0
ILMARTi |  70.80* |  73.20* |  77.57* | 153 | 48
Istella-S | NRGAM | 63.37 | 63.04 | 69.18 | 220 | 0
EBM | 57.75 | 59.67 | 66.90 | 220 | 0
EBMi | 62.29 | 63.05 | 69.71 | 220 | 50
ILMART | 66.99* | 66.29* | 72.36* | 106 | 0
ILMARTi |  69.06* |  68.22* |  74.04* | 106 | 50

Experimental results. Table 1 reports the results of the evaluation of ILMART,
ILMARTi and their state-of-the-art competitors. For all the datasets
considered, results show that ILMART achieves a statistically significant
improvement in terms of nDCG for almost all the cutoffs except the case of
nDCG@1 on Web30K. More importantly, ILMARTi always achieves the best
performance in terms of nDCG that always result in a statistically significant
improvement over NRGAM. The best improvements of ILMARTi over all the
competitors is obtained on Istella-S, where nDCG@10 improves of up to 6%
w.r.t. EBMi. On Web30K, we observe that ILMARTi improves nDCG@10 of up to 8%
w.r.t NRGAM, although the best competitor is, in this case, EBMi, whose
performance still lags behind ILMARTi. On Web30K, EBMi also achieves
statistically significant improvements in terms of nDCG with respect to NRGAM.
Even though EBMi is used as a pointwise ranker, these differences confirm the
efficacy of tree-based approaches in LtR and the need to include pairwise
interactions to learn more effective ranking functions.

Figure 2. nDCG by varying the number of interaction effects in ILMARTi on the
Web30K dataset. Figure 3. Plots of the four most important main and
interaction effects of ILMARTi on the Web30K dataset.

Furthermore, it is worth highlighting that our proposed techniques perform a
feature selection during the _main effects learning_ phase, while the
competitors leave the feature selection to the final user. This is an
important difference between ILMART and its competitors as a high number of
features used do not allow to easily understand the global behavior of the
model. For example, the Yahoo dataset comes with a large number of features,
i.e., $699$. ILMART and ILMARTi effectively reduce the initial set of $699$
features to $153$ (column $p$ in Table 1) with, at the same time, a
significant increase of the effectiveness of the learned model. In addition,
we point out that eventually ILMARTi may have a number of interaction effects
$K$ smaller than $50$ due to the training of the new $T$ trees which may not
include all the $K$ pairs selected.

We now investigate how the number of interactions added to ILMARTi can
influence its prediction accuracy. Fig. 2 compares the results of ILMARTi with
the ones of NRGAM by reporting their results in terms of nDCG. We observe that
the ranking accuracy of ILMARTi increases steeply when the first interaction
effects are added; then, the gain obtained by adding new interaction trees
saturates. This means that ILMARTi is able to significantly improve the
accuracy of an interpretable LtR model by adding a small number of interaction
effects, thus without hindering the overall interpretability of the model.
Even though the addition of interactions effects improves a lot the accuracy
of the model, we acknowledge that the difference from the full-complexity
model LambdaMART is still marked. For example, on Web30K an optimized black-
box LambdaMART model achieves a nDCG@10 close to 52%.

Finally, to show how the main and interaction effects can be explored, Fig. 3
presents three main effects and one interaction effects learned with ILMARTi
trained over Web30K. The three main effects selected are the ones associated
with the features with highest feature importance, i.e., PageRank (PR), Query-
url Click Count (QUCC), and Language Model for Information Retrieval (LIMIR),
while the interaction effect between URL Click Count (UCC) and QUCC is the one
with the highest average contribution for each value of $(x_{i},x_{j})$. The
2D plots and the heatmap are built by aggregating all the trees using the same
features, i.e., by representing the contribution of each $\tau_{j}(x_{j})$ and
each $\tau_{ij}(x_{i},x_{j})$. From the plots, the analyst can easily
understand the exact behavior of the model at varying values of each feature.
For example, taking into account the interaction effect, the higher the value
of UCC and QUCC the larger is the positive contribution of
$\tau_{ij}(x_{i},x_{j})$.

##  4\. Conclusions and Future Work

We presented ILMART, a novel and interpretable LtR solution based on
LambdaMART, that allows to train effective and intelligible models. We
experimentally assessed the performance of ILMART on three public datasets.
Results show that it achieves statistically significant improvements in terms
of nDCG of up to 8% w.r.t. the current state-of-the-art solution, i.e., NRGAM.
We also showed that ILMART is able to effectively exploit feature interactions
without hindering the overall interpretability of the model. With this work we
laid the foundations for the development of new accurate LtR models based on
boosting methods that are also interpretable.

As future work, we plan to investigate novel techniques based on boosting for
building effective and interpretable models for ranking and new approaches for
effectively identifying main and interaction effects during the learning of
the ranking models, with the main goal of filling the performance gap with the
black-box models.

## References

  * (1)
  * Burges et al. (2005) Chris Burges, Tal Shaked, Erin Renshaw, Ari Lazier, Matt Deeds, Nicole Hamilton, and Greg Hullender. 2005.  Learning to Rank Using Gradient Descent. In _Proceedings of the 22nd International Conference on Machine Learning_ (Bonn, Germany) _(ICML ’05)_. Association for Computing Machinery, New York, NY, USA, 89–96.  https://doi.org/10.1145/1102351.1102363
  * Castillo (2019) Carlos Castillo. 2019\.  Fairness and Transparency in Ranking.  _ACM SIGIR Forum_ 52, 2 (Jan. 2019), 64–71.  https://doi.org/10.1145/3308774.3308783
  * Chapelle and Chang (2011) Olivier Chapelle and Yi Chang. 2011.  Yahoo! Learning to Rank Challenge Overview. In _Proceedings of the Learning to Rank Challenge_ _(Proceedings of Machine Learning Research, Vol. 14)_ , Olivier Chapelle, Yi Chang, and Tie-Yan Liu (Eds.). PMLR, Haifa, Israel, 1–24.  https://proceedings.mlr.press/v14/chapelle11a.html
  * Cox (1984) D. R. Cox. 1984\.  Interaction.  _International Statistical Review / Revue Internationale de Statistique_ 52, 1 (1984), 1–24.  https://doi.org/10.2307/1403235 Publisher: [Wiley, International Statistical Institute (ISI)].
  * Fernando et al. (2019) Zeon Trevor Fernando, Jaspreet Singh, and Avishek Anand. 2019\.  A study on the Interpretability of Neural Retrieval Models using DeepSHAP. In _Proceedings of the 42nd International ACM SIGIR Conference on Research and Development in Information Retrieval_ _(SIGIR’19)_. Association for Computing Machinery, New York, NY, USA, 1005–1008.  https://doi.org/10.1145/3331184.3331312
  * Hastie and Tibshirani (1986) Trevor Hastie and Robert Tibshirani. 1986.  Generalized Additive Models.  _Statist. Sci._ 1, 3 (1986), 297 – 310.  https://doi.org/10.1214/ss/1177013604
  * Ke et al. (2017) Guolin Ke, Qi Meng, Thomas Finley, Taifeng Wang, Wei Chen, Weidong Ma, Qiwei Ye, and Tie-Yan Liu. 2017\.  LightGBM: A Highly Efficient Gradient Boosting Decision Tree. In _Proceedings of the 31st International Conference on Neural Information Processing Systems_ (Long Beach, California, USA) _(NIPS’17)_. Curran Associates Inc., Red Hook, NY, USA, 3149–3157.
  * Lou et al. (2012) Yin Lou, Rich Caruana, and Johannes Gehrke. 2012.  Intelligible models for classification and regression. In _Proceedings of the 18th ACM SIGKDD international conference on Knowledge discovery and data mining - KDD ’12_. ACM Press, Beijing, China, 150.  https://doi.org/10.1145/2339530.2339556
  * Lou et al. (2013) Yin Lou, Rich Caruana, Johannes Gehrke, and Giles Hooker. 2013\.  Accurate intelligible models with pairwise interactions. In _Proceedings of the 19th ACM SIGKDD international conference on Knowledge discovery and data mining - KDD ’13_. ACM Press, Chicago, Illinois, USA, 623.  https://doi.org/10.1145/2487575.2487579
  * Lucchese et al. (2017) Claudio Lucchese, Cristina Ioana Muntean, Franco Maria Nardini, Raffaele Perego, and Salvatore Trani. 2017.  RankEval: An Evaluation and Analysis Framework for Learning-to-Rank Solutions. In _Proceedings of the 40th International ACM SIGIR Conference on Research and Development in Information Retrieval_ (Shinjuku, Tokyo, Japan) _(SIGIR ’17)_. Association for Computing Machinery, New York, NY, USA, 1281–1284.  https://doi.org/10.1145/3077136.3084140
  * Lucchese et al. (2016) Claudio Lucchese, Franco Maria Nardini, Salvatore Orlando, Raffaele Perego, Fabrizio Silvestri, and Salvatore Trani. 2016.  Post-Learning Optimization of Tree Ensembles for Efficient Ranking. In _Proceedings of the 39th International ACM SIGIR Conference on Research and Development in Information Retrieval_ (Pisa, Italy) _(SIGIR ’16)_. Association for Computing Machinery, New York, NY, USA, 949–952.  https://doi.org/10.1145/2911451.2914763
  * Lundberg and Lee (2017) Scott M Lundberg and Su-In Lee. 2017.  A unified approach to interpreting model predictions.  _Advances in neural information processing systems_ 30 (2017).
  * Politou et al. (2018) Eugenia Politou, Efthimios Alepis, and Constantinos Patsakis. 2018\.  Forgetting personal data and revoking consent under the GDPR: Challenges and proposed solutions.  _Journal of Cybersecurity_ 4, 1 (03 2018).  https://doi.org/10.1093/cybsec/tyy001
  * Qin and Liu (2013) Tao Qin and Tie-Yan Liu. 2013.  Introducing LETOR 4.0 Datasets.  _CoRR_ abs/1306.2597 (2013).  http://arxiv.org/abs/1306.2597
  * Ribeiro et al. (2016) Marco Tulio Ribeiro, Sameer Singh, and Carlos Guestrin. 2016\.  ”Why Should I Trust You?”: Explaining the Predictions of Any Classifier. In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_. ACM, San Francisco California USA, 1135–1144.  https://doi.org/10.1145/2939672.2939778
  * Rudin (2019) Cynthia Rudin. 2019\.  Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead.  _Nature Machine Intelligence_ 1, 5 (May 2019), 206–215.  https://doi.org/10.1038/s42256-019-0048-x
  * Singh and Anand (2019) Jaspreet Singh and Avishek Anand. 2019.  EXS: Explainable Search Using Local Model Agnostic Interpretability. In _Proceedings of the Twelfth ACM International Conference on Web Search and Data Mining_. ACM, Melbourne VIC Australia, 770–773.  https://doi.org/10.1145/3289600.3290620
  * Smucker et al. (2007) Mark D. Smucker, James Allan, and Ben Carterette. 2007\.  A comparison of statistical significance tests for information retrieval evaluation. In _Proceedings of the sixteenth ACM conference on Conference on information and knowledge management - CIKM ’07_. ACM Press, Lisbon, Portugal, 623.  https://doi.org/10.1145/1321440.1321528
  * Verma and Ganguly (2019) Manisha Verma and Debasis Ganguly. 2019.  LIRME: Locally Interpretable Ranking Model Explanation. In _Proceedings of the 42nd International ACM SIGIR Conference on Research and Development in Information Retrieval_ _(SIGIR’19)_. Association for Computing Machinery, New York, NY, USA, 1281–1284.  https://doi.org/10.1145/3331184.3331377
  * Zhuang et al. (2021) Honglei Zhuang, Xuanhui Wang, Michael Bendersky, Alexander Grushetsky, Yonghui Wu, Petr Mitrichev, Ethan Sterling, Nathan Bell, Walker Ravina, and Hai Qian. 2021\.  Interpretable Ranking with Generalized Additive Models. In _Proceedings of the 14th ACM International Conference on Web Search and Data Mining_ _(WSDM ’21)_. Association for Computing Machinery, New York, NY, USA, 499–507.  https://doi.org/10.1145/3437963.3441796

  *[LtR]: Learning to Rank
  *[EBM]: Explainable Boosting Machine
  *[GAM]: Generalized Additive Models
  *[nDCG]: normalized Discounted Cumulative Gain
  *[PR]: PageRank
  *[QUCC]: Query-url Click Count
  *[LIMIR]: Language Model for Information Retrieval
  *[UCC]: URL Click Count