# Learning Interpretable Concept Groups in CNNs

Saurabh Varshneya1111Contact Author<EMAIL_ADDRESS>Antoine Ledent1
Robert A. Vandermeulen2 Yunwen Lei3 Matthias Enders4 Damian Borth5&Marius
Kloft1
1Technical University of Kaiserslautern, Germany
2Technical University of Berlin, Germany
3University of Birmingham, United Kingdom
4NPZ Innovation GmbH, Germany
5University of St.Gallen, Switzerland
{varshneya, ledent<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We propose a novel training methodology—Concept Group Learning (CGL)—that
encourages training of interpretable CNN filters by partitioning filters in
each layer into _concept groups_ , each of which is trained to learn a single
visual concept. We achieve this through a novel regularization strategy that
forces filters in the same group to be active in similar image regions for a
given layer. We additionally use a regularizer to encourage a sparse weighting
of the concept groups in each layer so that a few concept groups can have
greater importance than others. We quantitatively evaluate CGL’s model
interpretability using standard interpretability evaluation techniques and
find that our method increases interpretability scores in most cases.
Qualitatively we compare the image regions that are most active under filters
learned using CGL versus filters learned without CGL and find that CGL
activation regions more strongly concentrate around semantically relevant
features.††Code will be available at: https://github.com/srb-cv/cgl

##  1 Introduction

There is a great interest in understanding the hidden representations produced
by convolutional neural networks (CNNs). Understanding these representations
has significant implications theoretically (a better understanding of neural
networks in general) and practically (trustworthiness in AI). Toward this end,
many methods for interpreting hidden representations have been proposed. We
delineate two avenues of research, which include (1) visualizing the regions
of images that lead to high activation responses in a filter Yosinski et al.
(2014); Zhang et al. (2017); Mahendran and Vedaldi (2016) and (2) interpreting
the semantic meaning learned by a filter by finding its alignment to human-
interpretable visual concepts, such as object, scene, and color Bau et al.
(2017); Zhou et al. (2018); Zhang et al. (2020).

There exists much work that focuses on understanding the hidden
interpretations of CNNs, but there is a lack of work that actively influences
the learning process to favor interpretable representations.

To address this gap we propose a framework for CNNs that not only assesses
interpretability but induces it. Our approach is to induce a group structure
in CNNs at training, where groups of filters with a useful semantic meaning
emerge during the training process.

Figure 1: An illustration of our framework of learning concept groups in a
convolutional layer of a CNN.

We achieve this through carefully chosen regularization and auxiliary loss
functions. In our method, before any training is performed, we partition the
collection of filters in each layer into groups, which we term _concept
groups_. During training we promote filters within the same concept group to
learn similar features, so that, after training, they form a group of
correlated filters that jointly encode an abstract visual concept (e.g., an
object, scene, or color).

To learn these concept groups, we introduce three novelties to our training
objective (illustrated in Figure 1): 1\. We propose a new loss function
(entitled _group activation loss_) that encourages the filters within each
concept group to have highly overlapping activation regions for a given
sample. 2\. We propose another loss function (entitled _spatial loss_), which
promotes the activations of each filter to be concentrated around their mean
position. 3\. Using an $\ell_{2,1}$-block-norm regularizer on the concept
groups, we learn a sparse weighting of the concept groups. These weights can
be interpreted as the relative importance of each concept group. A key
property of our method is that it is completely unsupervised, automatically
finding concepts using zero side-information. In practice one would likely
include our losses in conjunction with some other task, e.g. classification,
to increase interpretability of that task, but none of the additional task
information, e.g. class labels, is used in our losses.

In addition to the training method described above, it can be observed that
interpretability of CNNs decreases when batch-normalization layers are
introduced Zhou et al. (2018). To remedy this we incorporate a variation of
batch-norm which works together with its regularizer and losses to achieve the
desired group structure in the hidden representations and enhance
interpretability.

We perform a quantitative analysis following the experimental setup by Bau et
al. (2017) and find that our method significantly improves interpretability.
We also analyze the filters of our model qualitatively by visualizing the
filters’ activations and find that our training setup yields representations
with noticeably more interpretable filters.

##  2 Related Work

We mention here a few existing works that are related to our training method
and outline the major differences with the existing approaches.

##### Group convolutions and structured sparsity.

Group convolutions were first used in Alexnet Krizhevsky et al. (2017) by
distributing filters into groups to train the model over multiple GPUs. It is
important to note that our concept group approach is different from group
convolutions, where there is no connection among the activation maps of
different partitioned groups. In contrast, concept groups’ activation maps are
concatenated at each layer before feeding them to the next layer (see Figure
3).

There exist methods that induce structured sparsity in CNNs Wen et al. (2016);
Li et al. (2017). Such methods show that inducing sparsity among groups of
filters can achieve state-of-the-art accuracy on various datasets with reduced
computational costs for CNNs. However, less is known about the effect of such
a regularization on the interpretability of networks.

##### Disentangled representations.

For improving hidden representations while training, a group of methods
endeavours to obtain _disentangled representations_. The aim of such methods
coincides with the main focus of our research. These methods aim to learn
better representations where each filter represents a clear and unique visual
concept. Zhang et al. (2018) introduce a training method which aims to train
each filter in the final convolutional layers of a CNN to represent an object
part. The activation loss from their method encourages a filter to activate
only for the training images that belong to the assigned category in the
training set.

Our own techniques also rely on applying a loss to the activation maps of
filters to inject priors into the CNN representation. However, there is a
major difference in how we apply the losses on activations in our study. In
Zhang et al. (2018) each filter is pushed towards a _known_ concept
corresponding to an object part. Ultimately the _concepts_ that this method
extracts are parts of the objects in the training classes. This is unlike our
method which can capture _arbitrary_ concepts (not only ones corresponding to
object parts) and without any supervision. Furthermore, the method from Zhang
et al. (2018) filters out activations of low magnitude entirely in the forward
pass. We argue that such hard constraints on the activation map can hamper the
representation power of the network. To retain representation power our method
applies soft constraints on the activations. This does not require any changes
in the forward pass.

The rest of this paper is organized as follows: in Section 3, we
mathematically define our algorithm and precise regularization procedure to
make the CNNs more interpretable. In Section 4.1, we explain the evaluation
techniques by Bau et al. (2017), which we later use to evaluate the
interpretability of the networks trained with our own training methods.
Finally, in Section 4.3, we compare the results of the experiments we
performed with the existing state-of-the-art methods.

##  3 Algorithm

Here we present our method which induces interpretable filters in all layers
of a CNN. We briefly outline our method here and describe them precisely in
the sequel. Central to our method is an initial grouping of each convolutional
layer’s filters into $G$ equal-size groups, e.g. for a collection of filters
in a layer indexed by $1,\ldots,F$ we divide the indices into $G$ groups
$(1,\ldots,f),(f+1,\ldots,2f),\ldots((G-1)f+1,\ldots,Gf)$ with $F=Gf$. During
training we aim to have a single group of filters (e.g. $(f+1,\ldots,2f)$)
correspond to a single concept. We assign multiple filters to the same concept
in the hopes of robustly capturing the concept in a way that would not be
possible with a single filter. The number of groups is a hyperparameter that
is determined before training. To make the method less restrictive, we can
leave some “free” filters outside the grouping. We consider each concept as a
block matrix and apply a general block norm over the collection of groups so
as to encourage structured sparsity Wen et al. (2016). The filters belonging
to each group are trained to correspond to similar visual concepts. To achieve
this we feed each filter’s activation map through a sigmoid to determine a
(soft) receptive field which represents the areas of the image where the
filter is active and then apply a loss to induce filters in the same concept
group to have receptive fields that largely overlap and are concentrated to a
connected region within an image (illustrated in Figure 2). We now precisely
describe the regularizers in the following.

Figure 2: The training setup for our auxiliary loss functions. At every layer
we obtain a soft receptive field by feeding a suitably standardized version of
the activation maps through a sigmoid. We apply constraints on the soft
receptive fields by our well-defined losses

###  3.1 Group Activation Loss $\operatorname{L_{g}}$

Here we describe our loss function that encourages filters in the same concept
group to learn similar concepts. To do this our loss function, which we denote
by $\operatorname{L_{g}}$, will encourage filters within the same concept
group to be active in roughly the same area for each input image. To determine
whether a filter is “active” in a region of the image, we first feed the
activation maps before applying any non-linearity (we call this the _pre-
activation map_) through a linear scaling (which is also learned and shared
between layers) followed by a sigmoid. We interpret the output of the sigmoid
as the probability of a filter being active at a location.

Let us denote the activation map of the $i^{\operatorname{th}}$ filter in
group $g$ in layer $l$ by $a_{ig}^{l}\in\mathbb{R}^{A_{l}}.$ Here, $A_{l}$
denotes the size of the pre-activation maps in layer $l$. A filter and its
activation maps are typically tensors but we simply consider a flattened
version of them for notational simplicity.

From our pre-activation map we utilize a sigmoid to yield a _soft-receptive
field_ , $\psi_{ig}^{l}\in\left(0,1\right)^{A_{l}}$, which is defined by the
following:

$\psi_{ig}^{l}=\sigma\left(P_{1}\frac{a_{ig}^{l}}{S_{ig}^{l}}+P_{2}\right),$
(1)

where $S_{ig}^{l}$ is the standard deviation of the $i^{\operatorname{th}}$
activation-map of group $g$, computed over the minibatch. Note that $P_{1}$
and $P_{2}$ are a pair of scaling parameters which can be learned during
training. These parameters are shared for all the activation maps in the
network. We now define the distance between two soft-receptive fields
$\psi_{1}$ and $\psi_{2}$ as

$\operatorname{\mathcal{D}}(\psi_{1},\psi_{2})=\frac{2\|\psi_{1}-\psi_{2}\|_{1}}{\|\psi_{1}\|_{1}+\|\psi_{2}\|_{1}+\|\psi_{2}-\psi_{1}\|_{1}},$
(2)

where $\|.\|_{1}$ denotes the $L^{1}$ norm: $\|x\|_{1}=\sum_{i}|x_{i}|$.

The motivation of this distance is to create a soft equivalent to one minus
the $\operatorname{IoU}$ (intersection over union). This score is frequently
used in computer vision to assess the performance of tracking or localization
algorithms Tychsen-Smith and Petersson (2018); Huang et al. (2019). To see the
equivalence, suppose that $\psi_{1}$ and $\psi_{2}$ are the indicator
functions of two sets $A$ and $B$. Then we have

$\displaystyle\mathcal{D}(\psi_{1},\psi_{2})=\frac{2\|\psi_{1}-\psi_{2}\|_{1}}{\|\psi_{1}\|_{1}+\|\psi_{2}\|_{1}+\|\psi_{1}-\psi_{2}\|_{1}}$
$\displaystyle=\frac{2\mu(A\Delta B)}{\mu(A)+\mu(B)+\mu(A\Delta
B)}=\frac{2\mu(A\Delta B)}{2\mu(A\cup B)}$ $\displaystyle=1-\frac{\mu(A\cap
B)}{\mu(A\cup B)}=1-\text{IOU},$ (3)

where we use the standard notation $A\Delta B:=(A\setminus B)\cup(B\setminus
A)=(A\cup B)\setminus(A\cap B)$ for the symmetric difference between $A$ and
$B$, and $\mu(.)$ denotes the Lebesgue measure or the cardinality in the case
of finite sets. The definition of $\psi$ must be altered slightly to work with
batch normalizaiton; we describe this in the following.

####  3.1.1 Definition of $\psi$ With Batch-normalization

If the network is trained with batch-norm, the pre-activations are
standardized to have mean $0$ and variance $1$. The input to the next layer is
then $\gamma_{ig}^{l}(a_{ig}^{l})+\beta_{ig}^{l}$, where $\beta_{ig}^{l}$ and
$\gamma_{ig}^{l}$ are the learned parameters of a batch-norm layer over the
$i^{\operatorname{th}}$ activation-map, of group $g$, as in standard batch
normalization Ioffe and Szegedy (2015). To modify the batch normalized values
to be compatible with our method we first define a variable $\tau_{ig}^{l}$ by

$\tau_{ig}^{l}=a_{ig}^{l}+\frac{(\beta_{ig}^{l})}{(\gamma_{ig}^{l})},$

which can be thought of as undoing batch normalization for the purposes of
finding a soft receptive field, since batch normalization is known to reduce
interpretability Bau et al. (2017). We can now feed $\tau_{ig}^{l}$ through a
sigmoid to obtain our soft receptive field as follows:

$\psi_{ig}^{l}=\sigma\left(P_{1}\tau_{ig}^{l}+P_{2}\right).$ (4)

####  3.1.2 Formula for $\operatorname{L_{g}}$

We can now define our group activation loss $\operatorname{L_{g}}$ to be the
following (note that the second term is new and is described in the next
paragraph)

$\displaystyle\frac{1}{r}\dfrac{\sum_{l,g}\sum_{(i,j)\in\mathcal{R}}\;\;\;\;\;\;\;\;\;\;\;\;\;2\|\psi_{{j}g}^{{l}}-\psi_{ig}^{l}\|_{1}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;}{\sum_{l,g}\sum_{(i,j)\in\mathcal{R}}\|\psi_{{j}g}^{{l}}\|_{1}+\|\psi_{ig}^{l}\|_{1}+\|\psi_{{j}g}^{{l}}-\psi_{ig}^{l}\|_{1}}+\cdots$
$\displaystyle\lambda_{2}\frac{1}{r}\dfrac{\sum_{l,g}\sum_{(i,j)\in\mathcal{R}}\;\;\;\;\;\;\;\;\;\;\;\;\;2\|\psi_{{j}g}^{{l}}-\psi_{ig}^{l+1}\|_{1}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;}{\sum_{l,g}\sum_{(i,j)\in\mathcal{R}}\|\psi_{{j}g}^{{l}}\|_{1}+\|\psi_{ig}^{l+1}\|_{1}+\|\psi_{{j}g}^{{l}}-\psi_{ig}^{l+1}\|_{1}}$
(5)

where $\mathcal{R}$ is a random subset of cardinality $r$ of
$\\{1,2,\ldots,N_{g}^{l}\\}\times\\{1,2,\ldots,N_{g}^{l}\\}$ that changes at
each iteration of the gradient descent procedure. $N^{l}_{g}$ denotes the
number of activation maps in group $g$ in layer $l$. We use this randomization
over $\mathcal{R}$ as an approximation of comparing all pairs of filters in
each concept group which is computationally expensive. Here we write
$\psi_{jg}^{l}$ for the vector composed of the components
$\psi_{jg}^{l}(x_{i})$ for samples $x_{i}$ in some minibatch. All the $L^{1}$
norms are computed over the minibatch (as well as the spatial dimension). The
value for $r$ can be set as a hyperparameter. We find that setting $r$ to
$3N_{g}^{l}$ works well in practice.

Each term in the summation that appears in $\operatorname{L_{g}}$ is a soft
version of one minus the intersection over union measure of similarity between
the receptive fields given by
$\operatorname{\mathcal{D}}(\psi_{ig}^{l},\psi_{jg}^{l})$. As some concepts
typically gradually change or die out in higher layers, we additionally
consider employing a sliding comparison, the second term in (5), where the
filters in group $g$ of layer $l$ are required to be similar to those in group
$g$ of layer $l+1$ but are not directly compared to those in group $g$ in
layers beyond $l+1$. We found that enforcing this similarity of concepts
throughout layers is too strong a requirement and simply causes the filters to
go to zero. Because of this we use $\lambda_{2}=0$ for all $l$ in our
experiments.

###  3.2 Spatial Loss $\operatorname{L_{s}}$

In this section we introduce a loss that aims to regularize each filter to
activate on a single, connected area of each image. Suppose for instance that
an image contains a car and a house both of which are relevant concepts to the
task to which the CNN is being applied. Since both concepts are abstract they
are likely to be represented near the final convolutional layers of a CNN
where the receptive fields are large. We would like to keep these concepts
separate between concept groups and prevent them from incorporating large
portions of the image corresponding to multiple abstract concepts. To achieve
this we enforce each filter to only be active on a small connected region on
images by penalizing activation responses that far from the center of
activation. The center of activation is computed over the soft receptive field
$\psi$ for each activation map. We define spatial loss
$\operatorname{L_{s}}(\psi)$ as the following:

$\operatorname{L_{s}}(\psi)=\frac{\sum_{j}\psi_{j}\|j-c(\psi)\|_{2}}{\sum_{j}\psi_{j}},\
\text{where }c(\psi)=\frac{\sum_{j}j\cdot\psi_{j}}{\sum_{j}\psi_{j}}$

and $j$ denotes the index of a activation of the soft receptive field ($\psi$)
of a filter. Here $c$ represents the mean position of activation of a soft
receptive field and $\operatorname{L_{s}}$ is analogous to the corresponding
variance. The loss is computed over all the activation maps of a filter
obtained in a minibatch.

Figure 3: The general block norm induces a group structure in the filters of a
layer, where few groups are dominantly learned over other groups. $\beta_{1}$,
$\beta_{2}$ and $\beta_{3}$ are the relevance factors of each group which are
learned while training a CNN

###  3.3 General Block Norm Regularizer $\operatorname{R_{bn}}$

In this section we define our general block norm regularizer
$\operatorname{R_{bn}}$ that replaces the conventional weight decay
regularizer. This changes the inductive bias of our training algorithm so that
it induces a group sparse structure over the concept groups. This regularizer
encourages some concept groups to be of greater importance than others,
thereby giving a notion of importance to concept groups and encouraging our
concept groups to be parsimonious. This allows us to avoid situations where,
for example, multiple concept groups capture the concept of “car.” Such sparse
models also yield reduced memory requirements and improved training speed as
side benefits Mocanu et al. (2018). An illustration the result of using this
regularization is shown in Figure 3. We will denote by $w_{gl}$ the matrix
containing all the weights of the filters of group $g$ in layer $l$. The
regularizer $\operatorname{R_{bn}}$ can then be written as :

$\operatorname{R_{bn}}=\sum_{l=1}^{L}\sum_{g=1}^{G}\|w_{gl}\|_{\operatorname{Fr}},$
(6)

where $L$ denotes the number of convolutional layers in the network, $G$ the
number of concept groups per convolutional layer, and
$\|\cdot\|_{\operatorname{Fr}}$ denotes the Frobenius norm i.e.,
$\|M\|_{\operatorname{Fr}}^{2}=\sum_{i,j}M_{i,j}^{2}$ for a matrix $M$.

An advantage of using the general block norm regularizer
$\operatorname{R_{bn}}$ is that we can obtain a vector of relevance factors
$\beta$ for our concept groups, as shown in Figure 3. These relevance scores
can be obtained by a summation of the norms as follows:

$\beta_{gl}=\|w_{gl}\|_{\operatorname{Fr}}\quad\text{and}\quad\beta_{l}=\sum_{g=1}^{G}\beta_{gl},$
(7)

where, $\beta_{gl}$ denotes a group’s relevance in layer $l$, and $\beta_{l}$
denotes the relevance of layer $l$ in a network.

##### Complete framework.

In summary our overall framework is obtained by combining the above mentioned
regularization strategy and the two loss functions. For a simplified view of
the combined method, we consider a CNN with weight parameters $w$ and we
denote the soft receptive fields obtained from the activation maps by $\psi$.
The final optimization target for an interpretable CNN can be defined as:

$\mathbf{L}=L_{d}(w)+\lambda_{bn}\operatorname{R_{bn}}(w)+\lambda_{g}\operatorname{L_{g}}(\psi)+\lambda_{s}\operatorname{L_{s}}(\psi),$
(8)

where, $L_{d}$ denotes the main loss for our task (e.g. cross entropy for
classification), $\operatorname{R_{bn}}$ is the general block norm
regularization, $\operatorname{L_{g}}$ and $\operatorname{L_{s}}$ are the
group activation loss and spatial loss respectively. $\lambda_{bn}$,
$\lambda_{g}$ and $\lambda_{s}$ are the weightings of the defined
regularization and auxiliary losses with respect to the main loss.

##  4 Experiments

In this section we explain the evaluation strategy we adopt for the
quantitative and qualitative analysis of the interpretability of CNNs, and
describe the results of our experiments performed on synthetic and real world
datasets.

###  4.1 Interpretability: Quantitative Assessment

In order to quantify the interpretability of a trained CNN, we exactly
replicate the evaluation on the Broden dataset as proposed by Bau et al.
(2017). The Broden dataset contains pixel-wise annotations of a broad range of
categories, which belong to one of six human-interpretable visual concepts. We
compute the alignment of a filter with a category by comparing its activation
maps for a set of images against the available ground truth using the same
threshold and scoring function described by Bau et al. (2017). The
interpretability of a layer can then be defined as the number of unique
categories assigned to the filters in that layer. Authors call this as the
number of unique detectors.

####  4.1.1 Evaluation on a Synthetic Dataset

In order to validate the efficacy of our methods in a controlled environment,
we constructed a synthetic dataset with obvious human-interpretable visual
concepts corresponding to shape and color: each image contains two randomly
positioned figures among the set {square, circle, triangle} with repetitions
allowed, and each shape is assigned a colour from the set {red, blue, green}.
Note that in principle, there are 15 high-level concepts in this paradigm: one
for each combination of colour and shape, and one for each colour and each
shape separately.

We consider two training settings: in the first case, we have one label for
each combination of shapes and colours (there are therefore 45 labels). In the
second setting, we simply assign the label one if a square is present and zero
otherwise (the problem is therefore binary). Therefore, in the first setting,
all the information contained in the high level concepts “shape, color, color-
shape” is contained in the label, whereas in the second situation, only part
of this information is encoded in the label.

For training, we use a simple CNN with $2$ convolutional layers, having $128$
and $256$ filters in the first and second layer respectively. We set the size
of filters to $3\times 3$ for both the layers. The evaluation is performed
using the same metrics as described above. We report the results of the
second, binary setting in Table 1. In the multiclass setting, the network is
able to recover all concepts very well regardless of regularization. However,
we find that in the binary setting, the addition of our regularizers has a
very strong positive effect on the recovery of the concepts: in the case where
we add all three regularizers, 14 out of 15 of all unique concepts are
recovered at the second layer, compared to only 10 when using weight decay.

This is a significant illustration of the fact that our method is able to
guide the network towards a more complete and disentangled representation of
all salient features of the data in the absence of low-level feature
annotations.

Regularizer | layer | No. of Unique Detectors
---|---|---
|  | color | shape | color-shape | Total
Weight-Decay | conv1 | 3 | 0 | 4 | 7
conv2 | 3 | 0 | 7 | 10
$\operatorname{R_{bn}}$ | conv1 | 3 | 2 | 4 | 9
conv2 | 3 | 2 | 7 | 12
$\operatorname{R_{bn}}+\operatorname{L_{g}}+\operatorname{L_{s}}$ | conv1 | 3 | 3 | 5 | 11
conv2 | 3 | 2 | 9 | 14
Table 1: Number of unique detectors in each layer while training on synthetic
data in a two class labels setup

####  4.1.2 Evaluation on the Real Datasets

We tested our methods on three well-known convolutional architectures named
Alexnet, Alexnet-B and VGG, where Alexnet-B is the version containing the
batch normalization layers added to the common Alexnet. All the networks are
trained from scratch on the two well-known image classification datasets
Places365 Zhou et al. (2017) and ImageNet Deng et al. (2009). The Places365
dataset is comparable to the ImageNet dataset, which contains 1.8 million
training images with 365 scene categories. For evaluation, we exactly
replicate the method from Bau et al. (2017), as explained above to compute the
interpretability of the trained models. For a baseline of interpretability, we
obtained the pretrained models for Imagenet and Places365 datasets from the
Pytorch and MIT CSAIL repositories respectively. Table 2 shows the comparison
of the interpretabilty scores for the models trained with the conventional
weight decay regularizer and our training method.

Dataset | Model | conv1 | conv2 | conv3 | conv4 | conv5 | score
---|---|---|---|---|---|---|---
Places365 | Alexnet | 6 | 12 | 22 | 36 | 54 | 130
Alexnet, ${(\operatorname{R_{bn}})}$ | 6 | 16 | 36 | 46 | 69 | 173
Alexnet, ${(\operatorname{R_{bn}}\operatorname{L_{g}}\operatorname{L_{s}})}$ | 8 | 17 | 35 | 50 | 75 | 185
Places365 | Alexnet-B | 6 | 13 | 27 | 43 | 59 | 148
Alexnet-B, ${(\operatorname{R_{bn}})}$ | 6 | 15 | 30 | 44 | 62 | 157
Alexnet-B, ${(\operatorname{R_{bn}}\operatorname{L_{g}}\operatorname{L_{s}})}$ | 5 | 17 | 30 | 51 | 63 | 166
Imagenet | Alexnet | 5 | 20 | 34 | 28 | 49 | 136
Alexnet, ${(\operatorname{R_{bn}})}$ | 5 | 16 | 39 | 32 | 45 | 137
Alexnet, ${(\operatorname{R_{bn}}\operatorname{L_{g}}\operatorname{L_{s}})}$ | 7 | 16 | 34 | 36 | 48 | 141
|  | conv3-3 | conv4-3 | conv5-1 | conv5-2 | conv5-3 |
Imagenet | VGG16 | 10 | 48 | 59 | 53 | 89 | 259
VGG16, ${(\operatorname{R_{bn}}\operatorname{L_{g}}\operatorname{L_{s}})}$ | 14 | 48 | 59 | 62 | 83 | 266
Places365 | VGG16 | 9 | 50 | 62 | 75 | 119 | 315
VGG16, ${(\operatorname{R_{bn}}\operatorname{L_{g}}\operatorname{L_{s}})}$ | 10 | 45 | 62 | 86 | 129 | 332

Table 2: The table above shows the sum of the number of unique detectors in
each layer of the network. We achieve a better interpretability across all the
tested networks and datasets.

####  4.1.3 Analysis over Concept-groups

For each predefined group, we compute a weighted sum over the number of
detectors found in the group and the $\operatorname{IoU}$ score. A group is
said to be aligned to a visual concept when the weighted sum is greater than a
threshold $T$, which was chosen by hand and is the same for all experiments.
One can obtain the relevance of such a concept-group from Equation 7. In
Figure 4, we show an analysis over the emergent concept-groups in the
convolutional layers of Alexnet.

Figure 4: Concept groups that emerge in different layers in a trained Alexnet.
The size of the circle represents the relevance of a concept group in a layer.

###  4.2 Interpretability: Qualitative Assessment

We use the method described in Zhou et al. (2014) for visualizing filters of
trained CNNs. To visualize the receptive field of a filter $f$, its activation
maps is computed over a set of images. The maximum activation (max-activation)
value for each activation map over all the pixels is extracted. The top $K$
images are then selected corresponding to the activation maps with $K$ largest
max-activation values. For each image, regions with “high” unit activations
are identified using a threshold $T_{f}$—taken from Bau et al. (2017). This
activated region is then scaled up to the image resolution and the filter
$f$’s receptive field can be visualized over the selected images. Overall,
this qualitative analysis gives a better understanding of a filter by focusing
on important areas in the image Zhou et al. (2014). We compare the
visualizations of the filters of a CNN trained with weight decay against a CNN
trained with our proposed losses in Figure 5.

Figure 5: The receptive fields of the filters in the last layer of Alexnet
trained with weight decay regularizer (in the top two rows) on Places365
dataset compared with the model trained with spatial loss
$\operatorname{L_{s}}$ (in the bottom two rows).

###  4.3 Accuracy Versus Interpretability

Most often, attempts to increase interpretability reduce the discrimination
power of a CNN, thereby compromising the validation accuracy of the trained
network. We show that by applying soft constraints as proposed in this paper,
the validation accuracy is only slightly reduced (less than 1 percent in each
case). Results in Table 2 show that most of the unique detectors emerge on the
last two convolutional layers of every network. Therefore we compare the ratio
of unique detectors (RUD), the number of the unique detectors divided by the
total number of filters in a layer, for the final two convolutional layers
along with the accuracy of the trained network in Table 3.

Dataset | Network | Top-1 acc. | Top-5 acc. | RUD
---|---|---|---|---
Places365 | Alexnet | 51.14% | 81.59% | 0.176
Alexnet, ${(\operatorname{R_{bn}}\operatorname{L_{g}}\operatorname{L_{s}})}$ | 50.46% | 80.60% | 0.244
| Alexnet-B* | 51.32% | 81.60% | 0.199
| Alexnet-B, ${(\operatorname{R_{bn}}\operatorname{L_{g}}\operatorname{L_{s}})}$ | 50.46% | 80.62% | 0.223
| VGG16 | 54.91% | 85.02% | 0.189
| VGG16, ${(\operatorname{R_{bn}}\operatorname{L_{g}}\operatorname{L_{s}})}$ | 54.91% | 85.20% | 0.210
Imagenet | Alexnet | 56.30% | 79.04% | 0.150
Alexnet, ${(\operatorname{R_{bn}}\operatorname{L_{g}}\operatorname{L_{s}})}$ | 55.09% | 77.93% | 0.164
| VGG16 | 71.79% | 90.45% | 0.139
| VGG16, ${(\operatorname{R_{bn}}\operatorname{L_{g}}\operatorname{L_{s}})}$ | 70.55% | 89.87% | 0.142
Table 3: A comparison of the Validation Accuracy and RUD Score on training
with different regularizers. The RUD score is the ratio of the total number of
unique detectors divided by the number of filters in the last two
convolutional layers of a CNN.

##  5 Conclusion

Our formulated losses and regularizer exploit the hidden structures of the
data via inducing a more prominent group structure among the filters of CNN.
Even simply replacing the conventional weight decay regularizer by our block
norm $\operatorname{R_{bn}}$, we can increase the interpretability of a
network without compromising with its discrimination power. We efficiently
compute the activation region of a filter over an image in the forward pass as
described in Section 3.1, and then constrain the activation regions so that
filters in the same group activate similar areas. This further induces a group
structure in all the layers of a CNN and enhances interpretability.

## Acknowledgements

The authors gratefully acknowledge support by the Carl-Zeiss Foundation, by
the German Research Foundation (DFG) award KL 2698/2-1, and by the Federal
Ministry of Science and Education (BMBF) awards 01IS18051A and 031B0770E.

## References

  * Bau et al. [2017] David Bau, Bolei Zhou, Aditya Khosla, Aude Oliva, and Antonio Torralba.  Network dissection: Quantifying interpretability of deep visual representations.  In Computer Vision and Pattern Recognition, 2017.
  * Deng et al. [2009] Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei.  Imagenet: A large-scale hierarchical image database.  In 2009 IEEE conference on computer vision and pattern recognition, pages 248–255. Ieee, 2009.
  * Huang et al. [2019] Yifeng Huang, Zhirong Tang, Dan Chen, Kaixiong Su, and Chengbin Chen.  Batching soft iou for training semantic segmentation networks.  IEEE Signal Processing Letters, 27:66–70, 2019.
  * Ioffe and Szegedy [2015] Sergey Ioffe and Christian Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  In Proceedings of the 32Nd International Conference on International Conference on Machine Learning - Volume 37, ICML’15, pages 448–456. JMLR.org, 2015.
  * Krizhevsky et al. [2017] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton.  Imagenet classification with deep convolutional neural networks.  Communications of the ACM, 60(6):84–90, 2017.
  * Li et al. [2017] Hao Li, Asim Kadav, Igor Durdanovic, Hanan Samet, and Hans Peter Graf.  Pruning filters for efficient convnets, 2017.
  * Mahendran and Vedaldi [2016] Aravindh Mahendran and Andrea Vedaldi.  Visualizing deep convolutional neural networks using natural pre-images.  International Journal of Computer Vision, 120(3):233–255, 2016\.
  * Mocanu et al. [2018] Decebal Constantin Mocanu, Elena Mocanu, Peter Stone, Phuong H Nguyen, Madeleine Gibescu, and Antonio Liotta.  Scalable training of artificial neural networks with adaptive sparse connectivity inspired by network science.  Nature communications, 9(1):1–12, 2018.
  * Tychsen-Smith and Petersson [2018] Lachlan Tychsen-Smith and Lars Petersson.  Improving object localization with fitness nms and bounded iou loss.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), June 2018.
  * Wen et al. [2016] Wei Wen, Chunpeng Wu, Yandan Wang, Yiran Chen, and Hai Li.  Learning structured sparsity in deep neural networks.  Advances in neural information processing systems, 29:2074–2082, 2016.
  * Yosinski et al. [2014] Jason Yosinski, Jeff Clune, Yoshua Bengio, and Hod Lipson.  How transferable are features in deep neural networks?  In Z. Ghahramani, M. Welling, C. Cortes, N. D. Lawrence, and K. Q. Weinberger, editors, Advances in Neural Information Processing Systems 27, pages 3320–3328. Curran Associates, Inc., 2014.
  * Zhang et al. [2017] Quanshi Zhang, Wenguan Wang, and Song-Chun Zhu.  Examining CNN representations with respect to dataset bias.  CoRR, abs/1710.10577, 2017.
  * Zhang et al. [2018] Quanshi Zhang, Ying Nian Wu, and Song-Chun Zhu.  Interpretable convolutional neural networks.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), June 2018.
  * Zhang et al. [2020] Quanshi Zhang, Xin Wang, Ying Nian Wu, Huilin Zhou, and Song-Chun Zhu.  Interpretable cnns for object classification, 2020.
  * Zhou et al. [2014] Bolei Zhou, Aditya Khosla, Àgata Lapedriza, Aude Oliva, and Antonio Torralba.  Object detectors emerge in deep scene cnns.  CoRR, abs/1412.6856, 2014.
  * Zhou et al. [2017] Bolei Zhou, Agata Lapedriza, Aditya Khosla, Aude Oliva, and Antonio Torralba.  Places: A 10 million image database for scene recognition.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 2017\.
  * Zhou et al. [2018] Bolei Zhou, David Bau, Aude Oliva, and Antonio Torralba.  Interpreting deep visual representations via network dissection.  IEEE transactions on pattern analysis and machine intelligence, 2018\.

  *[CNNs]: convolutional neural networks
  *[CNN]: convolutional neural network


# Two-Stage Learning of Highly Dynamic Motions
with Rigid and Articulated Soft Quadrupeds

Francecso Vezzi1,$\star$, Jiatao Ding1,$\star$,*, Antonin Raffin2, Jens
Kober1, Cosimo Della Santina1,2 This work is supported by the EU project
101016970 NI. 1Authors are with the Department of Cognitive Robotics, Delft
University of Technology, Building 34, Mekelweg 2, 2628 CD Delft, Netherlands
<EMAIL_ADDRESS><EMAIL_ADDRESS>C.DellaSantina@tudelft.nl). 2Authors
are with the Institute of Robotics and Mechatronics, German Aerospace Center
(DLR), 82234 Wessling, Germany (Antonin.Raffin@dlr.de). ⋆ These authors
contributed equally to this work. ∗ Jiatao Ding is the corresponding author.

###### Abstract

Controlled execution of dynamic motions in quadrupedal robots, especially
those with articulated soft bodies, presents a unique set of challenges that
traditional methods struggle to address efficiently. In this study, we tackle
these issues by relying on a simple yet effective two-stage learning framework
to generate dynamic motions for quadrupedal robots. First, a gradient-free
evolution strategy is employed to discover simply represented control
policies, eliminating the need for a predefined reference motion. Then, we
refine these policies using deep reinforcement learning. Our approach enables
the acquisition of complex motions like pronking and back-flipping,
effectively from scratch. Additionally, our method simplifies the
traditionally labour-intensive task of reward shaping, boosting the efficiency
of the learning process. Importantly, our framework proves particularly
effective for articulated soft quadrupeds, whose inherent compliance and
adaptability make them ideal for dynamic tasks but also introduce unique
control challenges.

##  I Introduction

Quadrupedal robots [1, 2, 3] promise to address real-world challenges due to
their potential advantages in terms of mobility and versatility. To improve
the performance and extend mission time [4], elastic and soft elements [5]
have been applied to quadrupedal robots in recent works, such as [6, 7, 8, 9].
Although stable locomotion can be achieved by off-the-shelf solvers,
accomplishing highly dynamic motion is still an open challenge due to the
under-actuated, hybrid, nonlinear dynamics.

Acrobatic motions have been achieved with quadrupedal robots by virtue of
trajectory optimization (TO) and model predictive control (MPC), including
jumping [10, 11, 12] and back-flip [13, 14]. However, these approaches are
computationally intensive and require accurate models. To guarantee the
solution when solving the TO problems, simplifications in problem formulations
are often needed, resulting in a conservative solution.

Reinforcement learning (RL) has emerged as a promising alternative to model-
based controllers for achieving dynamic motions. Although normal gaits like
trotting can be easily learned from scratch, learning acrobatic motions needs
extra setup, such as tedious reward shaping [15, 16, 17] and careful
curriculum design [18, 19, 20, 21], due to the highly nonlinear dynamics and
reward sparsity. To address this problem, the current RL frameworks assume the
availability of prior knowledge on how to solve the task in the form of a
reference motion/trajectory [22, 23, 24, 25, 26], a pre-existing controller
[27, 28] or both [29, 30]. This reliance on expert knowledge limits the
applicability of these methods as it may not be readily available and possibly
biases the learning towards sub-optimal policies.

Apart from motion generation for rigid robots, many works have argued that
introducing elasticity can yield better performance [31], which, however,
brings extra control changelings. Up to now, very limited work has been done
on generating highly dynamic quadrupedal motions exploiting parallel
elasticity. RL-based works such as [7] and [9] employ passive elasticity for
quadrupedal locomotion but do not report such acrobatic motions as done in
this work. The only exception is SpaceBok [8], which is, however, specifically
designed to work in a low-gravity environment.

To overcome these limitations, we propose here a two-stage strategy for
learning highly dynamic motions from scratch by combining RL with an evolution
strategy (ES), as illustrated in Fig. 1. To start with, an evolution stage is
employed to directly learn joint commands meeting the task demands, without
requiring an expert demonstration defined in advance. In this stage, we use
the augmented random search (ARS) [32] algorithm to learn a linear policy.
Then, we transfer the learned policy to a more expressive neural network via
deep reinforcement learning (DRL) [33]. We use proximal policy optimization
(PPO) [34] to refine the policy that is warm-started by mimicking the learned
motions in the first stage. Note that a similar two-stage method containing ES
and DRL has been proposed in [35] for quadrupedal locomotion. However, the ES,
specifically the natural evolution strategies algorithm [36], was used there
to search the parameters for a center pattern generator that generates
reference trajectories for the second-stage RL, and no dynamic motions like
ours were reported.

We validate111All the videos can be found at:
https://www.youtube.com/watch?v=liIeYU71a5w the approach by rich simulations
when learning versatile acrobatic motions, including jumping, pronking and
back-flip, see Fig. 2. We apply the proposed approach to both a rigid and an
articulated soft quadruped. We demonstrate that the proposed architecture
enables exploiting parallel elasticity in achieving highly dynamic motions.
Taking the jumping task as an example, the articulated soft quadrupedal robot
could jump 15.4% higher in the in-placed jumping task, and jump 23.1% further
in the forward jumping.

Figure 1: Two-stage learning procedure for the acrobatic motion with a flight
phase. Without defining reference motions, the (rigid and soft) quadrupedal
robot realizes (a) jumping in place and jumping forward, (b) pronking and (c)
back-flip. ‘ARS’ and ‘ARS+DRL’ separately represent the linear policies
generated by first-stage ES and the refined policy after the second-stage
retraining. Figure 2: Examples of highly dynamic motions that the robot could
learn using the proposed strategy. Red curves plot the CoM movements where the
arrows point to the movement directions.

##  II Methodology

We introduce basic notation in Sec. II-A, the two-stage procedure in a general
fashion in Sec. II-B. Acrobatic-specific details follow in Sec. II-C \-
yielding the general architecture shown in Fig. 1. Task-specific details will
be provided later in Sec. III.

###  II-A RL notation

The RL problem is formulated as a Markov decision process (MDP), described by
the tuple
$\mathcal{M}=\langle\mathcal{S},\mathcal{A},\mathcal{P},\mathcal{R},\gamma\rangle$
where $\mathcal{S}$ is the state space, $\mathcal{A}$ is the action space,
$\mathcal{P}$ is the transition probability, $\mathcal{R}$ is the reward
distribution, and $\gamma$ is the discount factor. The states
$(\bm{s}_{k}\in\mathcal{S})$, actions $(\bm{a}_{k}\in\mathcal{A})$, and
rewards $(r(\bm{s}_{k},\bm{a}_{k})\in\mathcal{R})$ in an episode of length $T$
constitute a rollout of the policy. In each rollout, the cumulative return is
$G_{T}=\sum_{k=1}^{T}\gamma^{k-1}r(\bm{s}_{k},\bm{a}_{k})$ . In the following,
we will assume that the whole state is observable, and thus we will refer to
elements in $\mathcal{S}$ as observations. The goal of RL is then to find an
approximation of the optimal policy ($\pi^{\star}$) such that the agent
achieves the maximal expected return, i.e.,
$\pi^{\star}=\arg\max\limits_{{\pi}}\mathbb{E}\left[G_{T}\mid\pi\right]$.

###  II-B Two-stage learning

####  II-B1 First step, learning a simple policy from scratch via evolution
strategy

Despite DRL’s success in sequential decision-making, it struggles with sparse-
reward tasks[37]. To address this, we use a gradient-free ES scheme,
specifically ARS, in the first stage to treat the RL tasks as black-box
optimization [38], where the linear policy is adopted222The work in [32]
argued that the linear policy is simple but enough to obtain competitive
results for RL tasks. A comparison with other representations on the jumping
task is found in Section IV-B.

$\small\pi_{\texttt{ars}}(\bm{s}_{k})=\mathbf{W}\,\bm{s}_{k}+\bm{b},$ (1)

where $\mathbf{W}:\mathcal{S}\rightarrow\mathcal{A}$ is a linear application,
and $\bm{b}\in\mathcal{A}$.

####  II-B2 Second step, refining and retraining via DRL

We then resort to the more expressive deep neural network (DNN) to refine the
policy. Using the DRL method, specifically PPO, we achieve the nonlinear
expression

$\small\pi_{\texttt{PPO}}(\bm{a}_{k}|\bm{s}_{k};{\bm{\theta}})=\mathrm{NN}(\bm{s}_{k};\bm{\theta}),$
(2)

where
$\mathrm{NN}:\mathcal{S}\times\mathbb{R}^{n}\rightarrow\Delta(\mathcal{A})$ is
a policy function that maps states to a probability distribution on the action
space, and $\bm{\theta}\in\mathbb{R}^{n}$ is the vector of the network
parameters.

At this stage, we start training a DNN by replicating the observation-action
pair provided by the first-stage policy, using the imitative reward function
$r^{\text{ini}}=w_{\mathrm{a}}\text{exp}^{(-w_{\mathrm{b}}||\bm{a}-\bm{a}^{\texttt{ars}}||^{2})}$,
with $\bm{a},\bm{a}^{\texttt{ars}}\in\mathcal{A}$ respectively being the
output of the DNN and the first-stage ARS given an observation, and
$w_{\mathrm{a}},w_{\mathrm{b}}\in\mathbb{R}$ being the hyperparameters.
Afterwards, we fine-tune the DNN using the task-specific reward ${}^{2}r$,
which we detail in section III.

###  II-C Learning acrobatic motions for quadrupeds

We define an acrobatic motion as a dynamic movement involving three (possibly
repeated) phases: (i) take-off, (ii) flight, (iii) landing. On top of what we
discussed in the previous section, we introduce a behavior selector mechanism.
This is a state machine that transitions through the three phases discussed
above, characterizing every acrobatic motion. The first phase is the take-off,
which uses (1). The second is the flight phase, which generates a constant
action equal to the homing pose in Fig. 3(b). The third phase passes the same
pose through a proportional-derivative (PD) controller thus reducing the
effective joint stiffness. This is sufficient to ensure a compliant enough
landing. The guarding conditions of the flight phase and landing phase
separately are $F_{\text{contact}}=0$, and $\dot{h}<0$ and
$F_{\text{contact}}>0$, where $F_{\text{contact}}$ is the total contact force
between the robot feet and the ground, and $\dot{h}$ is the vertical velocity.
We also add a low-pass filter to the action $a$ before feeding it into the
environment as policies learned with ARS can result in non-smooth control
signals. We remove these additions in the second stage of learning, relying on
the capability of the NN in (2) to learn the smooth actions.

Figure 3: Quadrupedal robot Go1 (a) and its PEA arrangement (b). In the homing
pose, we have initial height $z_{0}=0.32$m. We call $q_{1},q_{4}$ the calf
angles, $q_{2},q_{5}$ the thigh angles, and $q_{3},q_{6}$ the hip angles.

##  III Learning task formulation

This section defines the key elements of the RL tuple for learning versatile
acrobatic motions. We specifically focus on jumping, pronking, and back-flip
(see Fig. 2). Note that if not differently specified, the observation and
action spaces are the same during the two learning stages and the three tasks.

We conduct simulation experiments using the Go1 robot from UniTree and its
soft variant featuring parallel elastic actuation, as detailed in Appendix -A.
Both models are depicted in Fig. 3. Hence, when discussing specifics, we’ll
assume these are the robots in question, although our architecture is broadly
applicable to any quadrupedal system. Moreover, it is worth pointing out here
that we will focus on Sagittal motions - although the environment is fully 3D.

###  III-A Observation/state

These include the base (i.e., trunk) height ($h\in\mathbb{R}$), vertical
velocity ($\dot{h}\in\mathbb{R}$), and the pitch angle
($\theta\in\mathbb{R}$). Furthermore, the average333More precisely,
$q_{i}=(q_{i}^{\mathrm{left}}+q_{i}^{\mathrm{right}})/2$, where
$q_{i}^{\mathrm{left}},q_{i}^{\mathrm{right}}\in\mathbb{R}$ are the $i$-th
joint according to the numbering in Fig. 3 on the left and right sides of the
robot when cut by the Sagittal plane. joint angles ($\bm{q}\in\mathbb{R}^{6}$)
and velocities ($\dot{\bm{q}}\in\mathbb{R}^{6}$) are also added to the
observations. Finally, we include a Boolean flag $b\in\\{0,1\\}$ that
indicates whether the robot is in a flight phase or not. This observation is,
however, only available to $\pi_{\mathrm{PPO}}$. The resulting observations
vector is $\bm{s}_{k}=\\{h,\dot{h},\theta,\bm{q},\dot{\bm{q}},b\\}$. In
simulations, the white noise is added to the observations to emulate the
sensory noise.

###  III-B Action space

Assuming symmetry of the robot and of the task between the right and left
sides, we mirror the control action only learn the 3 joints in the front leg
and 3 joints in the rear leg.

We convert the policy outputs (i.e., $\bm{a}_{k}\in[-1,1]^{6}$) to joint angle
($\bm{q}^{\text{ref}}$) commands by a linear scaling defined such that the
maximum actions corresponds to the joint limits. which we feed into a PD
controller that generates motor torques as follows
$\tau_{i}=k_{\text{p}}(q^{\text{ref}}_{i}-q_{i})-k_{\text{d}}\dot{q}_{i},$
where $k_{\text{p}}\in\mathbb{R}$ and $k_{\text{d}}\in\mathbb{R}$ are the
proportional and derivative gains, respectively.

The pipeline of converting action to motor torque command is summarized by the
‘Action converter’ in Fig. 1.

###  III-C Episode design

In our settings, one episode is terminated when the time runs out or early
termination is triggered. In particular, the early termination is triggered if
one of the following is true

  * •

the robot falls, which we detect by checking the height, i.e., $h<0.1$ m,

  * •

the knees touch the ground, or links collide with each other,

  * •

the trunk is parallel enough to the ground, that we quantify by evaluating the
scalar product between the local $z$-axis and the global one and testing that
it is less than $0.85$.

The third condition is not present in the back-flip task.

###  III-D Rewards for the two stages

####  III-D1 Notation

We concisely refer to the exponential kernel as follows

$g_{l}^{t}:\mathbb{R}\to\mathbb{R},\quad g_{l}^{t}(x)\coloneqq
a_{l}^{t}e^{-b_{l}^{t}|x|}$

(3)

with $a_{l}^{t}$ and $b_{l}^{t}$ being hyperparameters. With this notation,
the kernel function for different tasks is encapsulated inside the subscript
and superscript, where the superscript ($t$) specifies the task, while the
subscript ($l$) refers to the input variable. Hereafter, we use $r$ to denote
$r(\bm{s}_{k},\bm{a}_{k})$ for brevity

####  III-D2 Reward shaping

To define the reward function for each task, we first normalize and clip the
maximal height ($\overline{h}\in\mathbb{R})$. For the jumping task, we also
normalize and clip the maximal jumping distance ($\overline{d}\in\mathbb{R}$).
That is

$\begin{aligned}
h_{n}=\text{clip}\left\\{0,{\overline{h}}/{h_{\max}},1\right\\},\quad\\!\\!d_{n}=\text{clip}\left\\{0,{\overline{d}}/{d_{\max}},1\right\\},\end{aligned}$

(4)

where $h_{\max}$ and $d_{\max}$ denote the maximal height and maximal jumping
distance that can be achieved by the robot, which we tuned heuristically.

We go now more in detail, taking the jumping task as an example. Pronking and
back-flip are similarly defined by making minor changes to the reward
functions. More details can be found in Appendix -B.

##### Reward function for the first stage

Since ES search can deal with the learning tasks with sparse rewards,
alleviating the efforts for rewarding shaping, we assume that the robot does
not get any reward during one episode until it ends. The other reason is that
dense reward might lead to unwanted behavior, whereas sparse reward makes it
simple to obtain the correct behavior. Thus, the reward for the first-stage
ARS training (${}^{1}r^{\text{jump}}$) is straightforwardly defined as

${}^{1}r^{\text{jump}}=\begin{cases}0&\text{when episode is running},\\\
r_{\text{bonus}}+r_{\text{end}}&\text{when episode terminates},\end{cases}$

(5)

where the ending reward $r_{\text{end}}$ is computed by evaluating the state
when the episode terminates, and the bonus reward $r_{\text{bonus}}$ is by
evaluating the jumping performance.

To achieve a jumping motion, we expect the robot to increase its normalized
height ($h_{n}$) as much as possible so as to take off from the ground.
Besides, to reduce the control efforts for posture adjustment in the air while
minimizing the landing recovery efforts after touchdown, we penalize the
maximal pitch angle ($\overline{\theta}\in\mathbb{R}$) the robot experiences
during one episode. Finally, one additional term is added to regulate the
normalized jumping distance ($d_{n}$). For jumping in place, this term
minimizes the forward distance while, for jumping forward, this term
encourages a large forward distance.

Reward for jumping in place are defined as

$\displaystyle\small r_{\text{end}}^{\text{jip}}$
$\displaystyle=\left[g_{{\theta}}^{\text{jip}}(\overline{\theta})+g_{d}^{\text{jip}}(d_{n})+c_{h}^{\text{jip}}\right]h_{n},$
(6) $\displaystyle r_{\text{bonus}}^{\text{jip}}$
$\displaystyle=\begin{cases}+\,b^{\text{jip}}h_{n}&\text{not early
termination},\\\
-\left[q^{\text{jip}}+m^{\text{jip}}h_{n}\right]&\text{otherwise},\end{cases}$

where $g_{\theta}^{\text{jip}}$ and $g_{d}^{\text{jip}}$ are computed using
(3). $c_{h}^{\text{jip}}$, $b^{\text{jip}}$, $q^{\text{jip}}$ and
$m^{\text{jip}}$ are the positive hyperparameters.

Reward for jumping forward are then defined as

$\displaystyle r_{\text{end}}^{\text{jf}}$
$\displaystyle=\left[g_{{\theta}}^{\text{jf}}(\overline{\theta})+c_{d}^{\text{jf}}d_{n}+c_{h}^{\text{jf}}\right]h_{n},$
$\displaystyle r_{\text{bonus}}^{\text{jf}}$
$\displaystyle=\begin{cases}+\,b^{\text{jf}}(h_{n}+d_{n})&\text{no early
termination},\\\
-\left[q^{\text{jf}}+m^{\text{jf}}(h_{n}+d_{n})\right]&\text{otherwise},\end{cases}$

where $g_{\theta}^{\text{jf}}$ is computed using (3). $c_{d}^{\text{jf}}$,
$c_{h}^{\text{jf}}$, $b^{\text{jf}}$, $q^{\text{jf}}$ and $m^{\text{jf}}$ are
the positive hyperparameters.

Note that the reward functions for jumping forward are quite similar to those
for jumping in place. The only difference is that the jumping distance appears
proportionally when computing $r_{\text{end}}^{\text{jf}}$ rather than inside
the kernel function such that a longer jumping distance is encouraged.
Besides, when computing the bonus for jumping forward
($r_{\text{bonus}}^{\text{jf}}$), we consider both the normalized maximal
jumping distance $d_{n}$ and the normalized maximal height $h_{n}$.

##### Reward function for the second stage

We discuss here the reward function used for refinement. The reader can refer
to Section II-B2 for the warm-starting reward $r^{\text{ini}}$. Then, the
dense refinement reward (${}^{2}r$) for retraining is

${}^{2}r^{\text{jump}}=r_{h}+r_{c}+r_{d}+r_{s}+r_{\theta}+{}^{2}r_{\text{bonus}}.$
(7)

In (7), the term $r_{h}$ is defined as

$r_{h}=\begin{cases}a_{h}h&\text{if }h_{\min}\leq h\leq h_{\max},\\\
0&\text{otherwise},\end{cases}$

(8)

with $a_{h}$ being a positive hyperparameter, which encourages the agent to
not exceed $h_{\max}$ when jumping and avoids rewarding it when not jumping at
all ($h<h_{\min}$). Note that here we use time-varying height $h$ instead of
the normalized maximal height $h_{n}$.

The term $r_{c}$ is used to penalize the average contact force
($f\in\mathbb{R}$) so as to encourage soft landing.

$r_{c}=\begin{cases}-a_{c}f&\text{if }f\geq f_{\min},\\\
0&\text{otherwise},\end{cases}$

(9)

where $f$ is computed by dividing total force ($F_{\text{contact}}$) by the
number of contact legs, $f_{\min}$ is the lower threshold, and $a_{c}$ is a
hyperparameter.

The term $r_{d}$ regulates jumping distance ($d\in\mathbb{R}$). Particularly,
for jumping in place, we minimize the jumping distance whereas, for jumping
forward, we maximize the jumping distance. To this end, we have

$\displaystyle r_{d}^{\text{jip}}=g_{d}^{\text{jip}}(d),\quad
r_{d}^{\text{jf}}=k_{d}^{\text{jf}}d.$ (10)

with $k_{d}^{\text{jf}}$ being a hyperparameter.

The term $r_{s}$ rewards a smooth torque profile
$r_{s}=g_{s}(\bm{\delta}_{\tau})$ with $\bm{\delta}_{\tau}\in\mathbb{R}^{6}$
being the differences between the torques applied at two consecutive time
steps.

The term $r_{\theta}$ minimizes the pitch angle, i.e.,
$r_{\theta}=g_{\theta}(\theta).$

Furthermore, similar to the first stage, a bonus is added at the end of one
episode. For both tasks, we define

$\begin{aligned} {}^{2}r_{\text{bonus}}^{\text{jip}}&=\begin{cases}0&\text{not
early termination},\\\ -m^{\text{jip}}\,\overline{h}&\text{otherwise}.\\\
\end{cases}\\\
{}^{2}r_{\text{bonus}}^{\text{jf}}&=\begin{cases}b^{\text{jf}}\,(\overline{h}+\overline{d})&\text{not
early termination},\\\ 0&\text{otherwise}.\\\ \end{cases}\\\ \end{aligned}$

(11)

with $m^{\text{jip}}$ and $b^{\text{jf}}$ being the hyperparameters.

##  IV Evaluation

The proposed architecture consistently succeeded in generating effective
acrobatic motions, as the one presented in Fig. 2 and in the video attachment.

In this section, we report an in-depth analysis of these performance,
including ablation studies and comparisons.

###  IV-A Technical details

We build the PyBullet environment [39] to emulate Go1’s motion. To mimic the
parallel spring in the training process, we add the resultant torque by spring
extensions (computed by the formula in Appendix -A) to each joint. To run the
experiments, we use the RL-Zoo training framework [40]. The neural network
implementing $\pi_{\mathrm{PPO}}$ is a two-layer DNN with 64 units each. We
use $\operatorname{tanh}$ as the activation function and set the learning rate
to be $2e^{-4}$. The batch size is 4096, and the discount factor $\gamma$ is
0.999. The clip range is 0.1.

Note: We will open-source all the code upon publication.

###  IV-B ES: linear policy vs nonlinear network

We first compare the linear policy with the NN-based nonlinear policies in ARS
training. Taking the in-placed jumping, for example, the learned curves using
different representations are plotted in Fig. 4. It demonstrates that the
linear policy (see ‘$\pi\\!\\!:\\!\\!\text{linear}$’ curve) is among the best
representations, obtaining the highest rewards. The comparable performance can
be achieved by the NN-based policy with one hidden layer of 64 units
(‘$\pi\\!\\!:\\!\\![64]$’ curve). Even though, the liner policy converges
faster with a smoother reward profile. Thus, it demonstrates that the linear
policy representation in the first-stage ARS training process is good enough
to achieve competitive results.

Figure 4: Reward profiles for the in-place jumping achieved by the ARS
algorithm when using different policy representations.

###  IV-C Learning explosive jumping motion

Using the hyperparameters in Appendix -C, the quadrupedal robot learns the
jumping motion. The first three columns in Fig. 5 present the results under
two stages, i.e., ‘ARS’ and ‘ARS+DRL’ (described in Fig. 1).

Figure 5: Jumping performance comparison. In the second column, a high reward
means a better performance, while a low contact force means a compliant
landing motion. Note the contact forces are divided by 1000N. In the third
column, the ‘height’ and ‘length’ separately denote the maximal jumping height
and jumping distance. In the fourth column, the learned CoM trajectories
(under ‘ARS+DRL’ policy) for a rigid robot and a soft robot are compared. ‘w/o
spring’ denotes the rigid case without springs engaged.

As can be seen from the first column of Fig. 5, both schemes can accomplish
two jumping tasks, i.e., jumping in place and jumping forward, successfully.
Among the two strategies, ‘ARS+DRL’ achieves the maximal mean reward (see the
green bar in the second column), accompanied by the larger jumping height for
both tasks (see the green bar in the third column). In addition, compared with
‘ARS’, a longer jumping distance is achieved by ‘ARS+DRL’ in the forward
jumping task, as can be seen from the yellow bar in the third column. Further
observations reveal that for in-placed jumping, ‘ARS+DRL’ achieves minimal
contact force, contributing to a compliant landing motion.

It should be mentioned that the first-stage ES itself (‘ARS’ in Fig. 5) could
accomplish two jumping tasks successfully, using the sparse reward defined in
Section III-D2. Considering the ‘ARS+DRL’ policy works better, in the
following section, we evaluate it by default.

TABLE I: Jumping performance with/without spring engaged | Jumping in place | Jumping forward
---|---|---
| spring | w/o spring | spring | w/o spring
Height [m] | 1.206$\pm$0.010 | 1.045$\pm$0.003 | 1.100$\pm$0.014 | 0.995$\pm$0.002
Distance [m] | 0.069$\pm$0.015 | 0.207$\pm$0.024 | 1.301$\pm$0.073 | 1.057$\pm$0.019
Force [N] | 799$\pm$189 | 871$\pm$228 | 796$\pm$211 | 800$\pm$369

###  IV-D Articulated soft quadruped vs rigid quadruped

Using the two-stage learning scheme, we train the control policies for both
the articulated soft robot (with parallel springs engaged) and the rigid robot
(without engaging parallel springs). To our knowledge, this is the first
report on learning quadrupedal jumping with parallel elasticity in a normal
gravity environment.

The fourth column of Fig. 5 compares center of mass (CoM) trajectories in both
cases (using ‘ARS+DRL’ policy), and Table I lists the meaning values of
jumping height, jumping distance and the largest landing force. It turns out
that a larger height (15.4% higher) and a longer distance (23.1% further) are
achieved when exploiting the parallel springs. Besides, parallel springs help
to gain smaller contact forces with fewer oscillations. More details can be
found in the attached video.

###  IV-E Beyond single jumping

Using the reward function defined in Appendix -B1 and Appendix -B2, the
proposed two-stage learning scheme enables learning the pronking and back-flip
motions, which are illustrated in Fig. 2. Note that no reference motion is
used here either. Fig. 6 presents the limit cycle behaviour when executing the
pronking motion. We find that the vertical CoM and pitch rotation would
converge to the limit cycle using the learned policy. Besides, we find that
without engaging parallel springs, the back-flip motion is hardly realized.

Figure 6: State evolutionary for pronking. Vertical motion (top) and pitch
rotation (bottom) converge to the limit cycle.

###  IV-F Environment randomization

To bridge the sim2real gap, we randomize the springs coefficient including
damping and stiffness, the total mass, robot mass distribution, CoM offset,
and the inertia tensor, as reported in Table II. In each episode, the
environment changes according to the random sampling.

To save space, we here only present the results of learning jumping motion.
Fig. 7 compares the jumping performance when using different strategies, where
‘ARS-r’ and ‘ARS+DRL-r’ separately denote the retrained policy using domain
randomization techniques [41]. For the in-placed jumping task, domain
randomization retraining helps to increase the success rate. In particular,
the two-stage retraining with dynamic randomization, i.e., ‘ARS+DRL-r’,
increases the jumping height. For forward jumping, the success rate is also
improved by domain randomization. One interesting thing is that the forward
jumping distance is decreased compared to the ‘ARS+DRL’ policy without domain
randomization. We guess it is because the learning scheme is prone to
guarantee the successful rate of forward jumping against environmental
uncertainties.

Figure 7: Jumping performances in randomized environments. ‘xx-r’ denotes the retrained policy with environment randomization. TABLE II: Environment randomization leg | mass | CoM | trunk | spring | spring
---|---|---|---|---|---
mass | offset | offset | mass | stiffness | damping
$\pm$20% | [0 4] [kg] |  $\pm$[0.2 0 0.2] [m] | accordingly |  $\pm$30% |  $\pm$30%

##  V Conclusions

In this work, we propose a two-stage approach for learning highly dynamic
quadrupedal motions from scratch. In the first stage, an evolution strategy is
used to generate the linear policy, without requiring defining a reference
trajectory in advance. Then, we refine and retrain the policy by deep
reinforcement learning, achieving higher performance. Compared with current
studies that focus on one specific task, it has been demonstrated that, by
slightly modifying the reward function, versatile acrobatic motions including
pronking and back-flip can be learned, with shared state space and observation
space.

Aside from the rigid case, we also apply the method to the articulated soft
robot with parallel elasticity. It turned out that the soft robot could still
learn highly dynamic motions using the proposed method. Furthermore,
exploitation of the parallel springs enhances the motion performance.

Future work will focus on hardware validation of these results on the real
robot.

###  -A Exploiting parallel elasiticity

We also explore parallel elasticity in achieving highly dynamic motions by
attaching parallel springs to the actuated joints of the Unitree Go1 robot
[3]. For the sake of simplicity, we here consider the springs attached to the
joints in the sagittal plane, i.e., thigh and calf joints (see Fig. 3). The
excessive torque on the $i$-th joint caused by the spring deformation
($\tau^{\text{s}}$) is computed as
$\tau^{\text{s}}_{i}=k_{i}\left(q_{i}^{\text{ref}}-q^{\text{0}}_{i}\right)+c_{i}\dot{q}_{i}^{\text{ref}}$
with $k_{i}\in\mathbb{R}$ and $c_{i}\in\mathbb{R}$ separately being the spring
stiffness and damping, $q_{i}^{\text{ref}}\in\mathbb{R}$,
$\dot{q}_{i}^{\text{ref}}\in\mathbb{R}$, and $q^{\text{0}}_{i}\in\mathbb{R}$
separately denoting the commanded angle, angular velocity) and rest angle of
the $i$-th joint.

###  -B Rewards for pronking and back-flip

####  -B1 Reward for pronking

To begin with, we evaluate single jump by

$p_{\text{sj}}=w_{h}h_{n}+w_{d}d_{n}.$ (12)

Then, we define the performance array for $n$ jumps during one episode as
$p_{j}=\left[p_{\text{sj}}^{1},p_{\text{sj}}^{2},\dots,p_{\text{sj}}^{n}\right]$.
Given $p_{j}$, we can compute the average performance ${p}_{j}^{\text{avg}}$
and the performance entropy $S(p_{j})$ (considering
$p_{\text{sj}}^{i}\in\left[0,1\right]$ and
$\sum_{i=1}^{n}p_{\text{sj}}^{i}=1$). The sparse reward used for the first-
stage learning is then defined as

$\begin{aligned}
&r_{\text{avg}}={p}_{j}^{\text{avg}}\left\\{g_{\theta}(\overline{\theta})+w_{t}t_{n}+w_{s}S(p_{j})+k\right\\},\\\
&{}^{1}r^{\text{pronking}}=w_{\text{avg}}r_{\text{avg}}+w_{\text{max}}p_{\text{sj}}^{\max}+w_{c}c_{j}+b,\end{aligned}$

(13)

where $c_{j}$ is the number of jumps whose performance is above a certain
threshold, $b$ is a bonus for a successful termination, $t_{n}$ is the ratio
between the time the episode ends and the maximal episode time, $k$ is a
constant and $p_{\text{sj}}^{\max}$ is the maximal element of the array
$p_{j}$. $w_{\times}$ denotes the hyperparameter.

For the second-stage training, the reward function is similar to the one
already defined in (7), added by two additional terms. One term ($r_{p}$) is
used to minimize the actuators’ energy consumption, and the other term
($r_{j}$) provides a reward each time the robot successfully performs a jump,
which are separately defined as

$\begin{aligned}
&r_{p}=g_{p}\left(E\right),\quad&r_{j}=w_{j}{p}_{j}^{\text{avg}}\left[g_{s}\left(S\left(p_{j}\right)\right)+l_{j}\right],\end{aligned}$

(14)

where $E$ is the energy cost in one episode, $l_{j}$ is a constant reward, and
the weight $w_{j}$ is $0$ if the jump performance is below a certain
threshold.

####  -B2 Reward for back-flipping

Denoting $\theta_{n}\in\mathbb{R}$ as the normalized maximum pitch angle in
one episode (computed by (4)), the reward function for the first-stage ARS
training is defined as

${}^{1}r^{\text{flip}}=w_{h}\overline{h}+w_{\theta}\theta_{n}+w_{\text{h}\theta}(\overline{h}\theta_{n})+b,$
(15)

where $b$ is a bonus if the episode ends successfully, $w_{\times}$ denotes
the hyperparameter.

The reward for the second-stage DRL training is similar to the one defined in
(7). One exception is the reward with regard to the jumping distance ($r_{d}$
in (7)) is discarded and the one with regard to the pitch angle is
proportional to it, i.e., $r_{\theta}=g_{\theta}(\theta)$. Another point is
that we change $r_{\text{bonus}}$ in (7) to be

$r_{\text{bonus}}=\begin{cases}w_{\text{bonus}}(\overline{h}\theta_{n})&\text{not
early termination},\\\ 0&\text{otherwise}.\end{cases}$

(16)

###  -C Hyperparameters for learning jumping motion

Hyperparameters for the first-stage and second-stage training are separately
listed in Table III and Table IV.

TABLE III: Reward penalties for the first-stage ARS training. Jumping | $h_{f}^{\text{jip}}$ | $a_{\theta}^{\text{jip}}$ | $b_{\theta}^{\text{jip}}$ | $a_{d}^{\text{jip}}$ | $b_{d}^{\text{jip}}$ | $c_{h}^{\text{jip}}$ | $b^{\text{jip}}$ | $q^{\text{jip}}$ | $m^{\text{jip}}$
---|---|---|---|---|---|---|---|---|---
in place | 0.9 | 0.3 | 0.0225 | 0.05 | 0.05 | 0.7 | 0.1 | 0.08 | 0.064
Jumping | $h_{f}^{\text{jf}}$ | $a_{\theta}^{\text{jf}}$ | $b_{\theta}^{\text{jf}}$ | $d_{f}^{\text{jf}}$ | $c_{d}^{\text{jf}}$ | $c_{h}^{\text{jf}}$ | $b^{\text{jf}}$ | $q^{\text{jf}}$ | $m^{\text{jf}}$
forward | 0.3 | 0.25 | 0.0225 | 1.3 | 0.5 | 25 | 0.1 | 0.08 | 0.096
TABLE IV: Reward penalties for the second-stage PPO training. Jumping | $a_{h}^{\text{jip}}$ | $h_{\min}^{\text{jip}}$ | $h_{\max}^{\text{jip}}$ | $a_{c}^{\text{jip}}$ | $f_{\min}^{\text{jip}}$ | $a_{d}^{\text{jip}}$
---|---|---|---|---|---|---
0.01 | 0.29 | 1.3 | $1.5e^{-4}$ | 800 | $6.5e^{-4}$
in place | $b_{d}^{\text{jip}}$ | $a_{s}^{\text{jip}}$ | $b_{s}^{\text{jip}}$ | $a_{\theta}^{\text{jip}}$ | $b_{\theta}^{\text{jip}}$ | $m^{\text{jip}}$
| 40 | $3e^{-3}$ | 0.1 | $4.2e^{-3}$ | 26 | 0.25
Jumping | $a_{h}^{\text{jf}}$ | $h_{\min}^{\text{jf}}$ | $h_{\max}^{\text{jf}}$ | $a_{c}^{\text{jf}}$ | $f_{\min}^{\text{jf}}$ | $k_{d}^{\text{jf}}$
$6.5e^{-3}$ | 0.29 | 1.1 | $1.2e^{-4}$ | 800 | $1.52e^{-2}$
forward | $d_{\max}^{\text{jf}}$ | $a_{s}^{\text{jf}}$ | $b_{s}^{\text{jf}}$ | $a_{\theta}^{\text{jf}}$ | $b_{\theta}^{\text{jf}}$ | $b^{\text{jf}}$
| 1.3 | $3e^{-3}$ | 0.1 | $4.2e^{-3}$ | 26 | 0.025

## References

  * [1] M. Hutter, C. Gehring, D. Jud, A. Lauber, C. D. Bellicoso, V. Tsounis, J. Hwangbo, K. Bodie, P. Fankhauser, M. Bloesch et al., ”Anymal-a highly mobile and dynamic quadrupedal robot,” in _Proc. IEEE/RSJ Int. Conf. Intell. Robot. Syst._ , 2016, pp. 38–44.
  * [2] B. Katz, J. Di Carlo, and S. Kim, “Mini cheetah: A platform for pushing the limits of dynamic quadruped control,” in _Proc. IEEE Int. Conf. Robot. Autom._ , 2019, pp. 6295–6301.
  * [3] U. Robotics, “Unitree robot go1: https://m.unitree.com/products/go1,” 2021\.
  * [4] F. Angelini, P. Angelini, C. Angiolini, S. Bagella, F. Bonomo, M. Caccianiga, C. Della Santina, D. Gigante, M. Hutter, T. Nanayakkara et al., “Robotic monitoring of habitats: the natural intelligence approach,” _IEEE Access_ , 2023.
  * [5] C. Della Santina, M. G. Catalano, A. Bicchi, M. Ang, O. Khatib, and B. Siciliano, “Soft robots,” _Ency. Robot._ , vol. 489, 2021.
  * [6] A. Badri-Sprowitz, A. Aghamaleki Sarvestani, M. Sitti, and M. A. ¨ Daley, “Birdbot achieves energy-efficient gait with minimal control using avian-inspired leg clutching,” _Sci. Robot._ , vol. 7, no. 64, p.eabg4055, 2022.
  * [7] F. Bjelonic, J. Lee, P. Arm, D. Sako, D. Tateo, J. Peters, and M. Hutter, “Learning-based design and control for quadrupedal robots with parallel-elastic actuators,” _IEEE Robot. Autom. Lett._ , vol. 8, no. 3, pp. 1611–1618, 2023.
  * [8] P. Arm et al., “Spacebok: A dynamic legged robot for space exploration,” in _Proc. IEEE Int. Conf. Robot. Autom._ , 2019, pp. 6288–6294.
  * [9] A. Raffin, D. Seidel, J. Kober, A. Albu-Schaffer, J. Silv ¨ erio, and F. Stulp, “Learning to exploit elastic actuators for quadruped locomotion,” _arXiv preprint arXiv:2209.07171_ , 2022.
  * [10] Q. Nguyen, M. J. Powell, B. Katz, J. Di Carlo, and S. Kim, “Optimized jumping on the mit cheetah 3 robot,” in _Proc. IEEE Int. Conf. Robot. Autom._ , 2019, pp. 7448–7454.
  * [11] C. Nguyen, L. Bao, and Q. Nguyen, “Continuous jumping for legged robots on stepping stones via trajectory optimization and model predictive control,” in _Proc. IEEE Conf. Decis. Control_ , 2022, pp. 93–99.
  * [12] Z. Song, L. Yue, G. Sun, Y. Ling, H. Wei, L. Gui, and Y.-H. Liu, “An optimal motion planning framework for quadruped jumping,” in _Proc. IEEE/RSJ Int. Conf. Intell. Robot. Syst._ , 2022, pp. 11 366–11 373.
  * [13] Y. Ding, A. Pandala, C. Li, Y.-H. Shin, and H.-W. Park, “Representation-free model predictive control for dynamic motions in quadrupeds,” _IEEE Trans. Robot._ , vol. 37, no. 4, pp. 1154–1171, 2021.
  * [14] G. Garc´ıa, R. Griffin, and J. Pratt, “Time-varying model predictive control for highly dynamic motions of quadrupedal robots,” in _Proc. IEEE Int. Conf. Robot. Autom._ , 2021, pp. 7344–7349.
  * [15] J. Hwangbo, J. Lee, A. Dosovitskiy, D. Bellicoso, V. Tsounis, V. Koltun, and M. Hutter, “Learning agile and dynamic motor skills for legged robots,” _Sci. Robot._ , vol. 4, no. 26, p. eaau5872, 2019.
  * [16] Z. Fu, A. Kumar, J. Malik, and D. Pathak, “Minimizing energy consumption leads to the emergence of gaits in legged robots,” _arXiv preprint arXiv:2111.01674_ , 2021.
  * [17] 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.
  * [18] [18] Z. Tang, D. Kim, and S. Ha, “Learning agile motor skills on quadrupedal robots using curriculum learning,” in _Proc. Int. Conf. Robot. Intell. Tech. Appl._ , vol. 3, 2021.
  * [19] [19] G. B. Margolis, G. Yang, K. Paigwar, T. Chen, and P. Agrawal, “Rapid locomotion via reinforcement learning,” _arXiv preprint arXiv:2205.02824_ , 2022.
  * [20] [20] A. Iscen, G. Yu, A. Escontrela, D. Jain, J. Tan, and K. Caluwaerts, “Learning agile locomotion skills with a mentor,” in _Proc. IEEE Int. Conf. Robot. Autom._ , 2021, pp. 2019–2025.
  * [21] S. Chen, B. Zhang, M. W. Mueller, A. Rai, and K. Sreenath, “Learning torque control for quadrupedal locomotion,” _arXiv preprint arXiv:2203.05194_ , 2022.
  * [22] G. B. Margolis, T. Chen, K. Paigwar, X. Fu, D. Kim, S. Kim, and P. Agrawal, “Learning to jump from pixels,” _arXiv preprint arXiv:2110.15344_ , 2021.
  * [23] Y. Fuchioka, Z. Xie, and M. Van de Panne, “Opt-mimic: Imitation of optimized trajectories for dynamic quadruped behaviors,” in _Proc. IEEE Int. Conf. Robot. Autom._ , 2023, pp. 5092–5098.
  * [24] M. Bogdanovic, M. Khadiv, and L. Righetti, “Model-free reinforcement learning for robust locomotion using demonstrations from trajectory optimization,” _arXiv preprint arXiv:2107.06629_ , 2021.
  * [25] C. Li, S. Blaes, P. Kolev, M. Vlastelica, J. Frey, and G. Martius, “Versatile skill control via self-supervised adversarial imitation of unlabeled mixed motions,” _arXiv preprint arXiv:2209.07899_ , 2022.
  * [26] C. Li, M. Vlastelica, S. Blaes, J. Frey, F. Grimminger, and G. Martius, “Learning agile skills via adversarial imitation of rough partial demonstrations,” in _Conference on Robot Learning. PMLR_ , 2023, pp. 342–352.
  * [27] L. Smith, J. C. Kew, T. Li, L. Luu, X. B. Peng, S. Ha, J. Tan, and S. Levine, “Learning and adapting agile locomotion skills by transferring experience,” _arXiv preprint arXiv:2304.09834_ , 2023.
  * [28] Y. Yang, X. Meng, W. Yu, T. Zhang, J. Tan, and B. Boots, “Continuous versatile jumping using learned action residuals,” in _Learning for Dynamics and Control Conference. PMLR_ , 2023, pp. 770–782.
  * [29] Y. Yang, G. Shi, X. Meng, W. Yu, T. Zhang, J. Tan, and B. Boots, “Cajun: Continuous adaptive jumping using a learned centroidal controller,” _arXiv preprint arXiv:2306.09557_ , 2023.
  * [30] Z. Li, X. B. Peng, P. Abbeel, S. Levine, G. Berseth, and K. Sreenath, “Robust and versatile bipedal jumping control through multi-task reinforcement learning,” _arXiv preprint arXiv:2302.09450_ , 2023.
  * [31] N. Kashiri, A. Abate, S. J. Abram, A. Albu-Schaffer, P. J. Clary, M. Daley, S. Faraji, R. Furnemont, M. Garabini, H. Geyer et al., “An overview on principles for energy efficient robot locomotion,” _Front. Robot. AI_ , vol. 5, p. 129, 2018.
  * [32] H. Mania, A. Guy, and B. Recht, “Simple random search of static linear policies is competitive for reinforcement learning,” _Adv. Neural Inf. Process._ , vol. 31, 2018.
  * [33] V. Franc¸ois-Lavet, P. Henderson, R. Islam, M. G. Bellemare, J. Pineau et al., “An introduction to deep reinforcement learning,” _Found. Trends Mach. Learn._ , vol. 11, no. 3-4, pp. 219–354, 2018.
  * [34] J. Schulman, F. Wolski, P. Dhariwal, A. Radford, and O. Klimov, “Proximal policy optimization algorithms,” _arXiv preprint arXiv:1707.06347_ , 2017.
  * [35] H. Shi, B. Zhou, H. Zeng, F. Wang, Y. Dong, J. Li, K. Wang, H. Tian, and M. Q.-H. Meng, “Reinforcement learning with evolutionary trajectory generator: A general approach for quadrupedal locomotion,” _IEEE Robot. Autom. Lett._ , vol. 7, no. 2, pp. 3085–3092, 2022.
  * [36] T. Salimans, J. Ho, X. Chen, S. Sidor, and I. Sutskever, “Evolution strategies as a scalable alternative to reinforcement learning,” _arXiv preprint arXiv:1703.03864_ , 2017.
  * [37] A. Y. Majid, S. Saaybi, V. Francois-Lavet, R. V. Prasad, and C. Verhoeven, “Deep reinforcement learning versus evolution strategies: a comparative survey,” _IEEE Trans. Neural Netw. Learn._ , 2023.
  * [38] N. Hansen, D. V. Arnold, and A. Auger, “Evolution strategies,” _Springer handbook of computational intelligence_ , pp. 871–898, 2015.
  * [39] E. Coumans and Y. Bai, “Pybullet, a python module for physics simulation for games, robotics and machine learning,” _http://pybullet.org_ , 2016–2019.
  * [40] A. Raffin, “Rl baselines3 zoo,” https://github.com/DLR-RM/ rl-baselines3-zoo, 2020.
  * [41] Z. Xie, X. Da, M. Van de Panne, B. Babich, and A. Garg, “Dynamics randomization revisited: A case study for quadrupedal locomotion,” in _Proc. IEEE Int. Conf. Robot. Autom._ , 2021, pp. 4955–4961.


# Validation of x-ray line-profile analysis of extended defects in deformed
crystals

C P Race<EMAIL_ADDRESS>T Ungar  Department of Materials,
University of Manchester, Manchester, M13 9PL, United Kingdom  G Ribárik
Department of Materials Physics, Eötvös Loránd University Budapest, H-1117
Pázmány P. sétány 1/A, Hungary

###### Abstract

Quantitative measurements of extended defects in crystalline materials are
important in understanding material behaviour. X-ray line profile analysis
provides a complement to direct counting in the electron microscope, but is an
indirect method and requires validation. Previous studies have focused on
comparing x-ray analysis to electron microscopy results. Instead, we use
simulated defective material with known defect content and apply line profile
analysis to calculated diffraction profiles to directly show that line profile
analysis can reliably quantify dislocations and stacking faults.

Extended defects in crystalline materials, such as dislocations and stacking
faults, strongly influence material properties and play a key role in many
failure mechanisms. To properly understand material behaviour, we therefore
need reliable ways of quantifying and characterising these defects. One common
approach is to directly image defects in the transmission electron microscope
and count them. Though highly successful, this approach has several
disadvantages: it is time consuming and expensive to prepare the thin foils
needed for analysis; this preparation can strongly affect the defects that we
wish to measure; not all the defects will be large enough to see; high defect
densities make individual defects difficult to identify; and the volume of
material analysed in a typical foil is only very small.

An alternative approach is to use x-ray _line profile analysis_. In this
approach, we examine the one-dimensional profile of diffracted intensity as a
function of angle, arrived at via an integration around the Debye-Scherrer
rings formed by diffraction of a highly mono-chromated beam of x-rays. Such
profiles are commonly generated from powder or polycrystalline samples, but
with a bright enough source of x-rays (e.g. from a synchrotron) profiles from
individual grains can be measured. The _positions_ of the peaks in the line
profiles correspond to the _spacing_ of sets of planes in the sample and so
can be used to determine crystal structure and macro strains. Line profile
analysis takes this a step further by making use of the fact that the _shapes_
of the peaks are affected by the _distortion_ of the corresponding set of
planes and so contain information about the types and numbers of defects
present.

Because line profile analysis requires only simple sample preparation and can
be used to quickly analyse a relatively large volume of material, it provides
an important complement to direct counting of defects in the electron
microscope. We need to remember, however, that it is an indirect method of
quantifying defect content. It begins with physical models of the distortions
created by individual defects and considers the aggregate effect that a
population of such defects will have on the shapes of the peaks in a line
profile. This model line profile is thus a function of parameters representing
the nature, density and arrangement of defects within the material and the
technique proceeds by optimising the values of these parameters to achieve a
model profile that best fits the experimental profile. We are thus attempting
to infer details of complex features of the sample from a relatively simple
one-dimensional signal and the task would be hopeless, were it not the case
that the response of the peak shapes varies from peak to peak, with
diffraction order and with distance from the peak centre, in different ways
for different defects and for different types and arrangements of
dislocations.

The physical basis for line profile analysis is sound. The distortion due to a
dislocation, for example, is well described by theory at various levels of
physics, including for the case of a general burgers vector and line
direction, and for anisotropic elasticity Bacon:1980fk . However, given the
complexity of the process, an empirical validation is required. Previous
attempts to validate the results of line profile analysis have relied on
comparisons with counts in the electron microscope Ungar:1984aa , but this is
a somewhat indirect approach: the electron microscope does not give the full
picture. Here we undertake an alternative, more direct, validation.

Rather than consider an experimental sample of material, we instead take a
sample of simulated material to which we can apply line profile analysis and
for which we know the precise defect content. We proceed by taking an
atomistic simulation cell containing $\sim 4$ million atoms in a face-centred
cubic arrangement. We use an embedded atom method model for copper
Ackland:1987hs  and randomly add six $1/2\langle 110\rangle$-type dislocation
loops of $12\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ radius to act as “seeds”
for plastic deformation. Using the classical molecular dynamics software
Lammps Plimpton:1995fv  we briefly anneal the simulation cell for
$10\,\mathrm{ps}$ at $300\,\mathrm{K}$ then apply a uniaxial distortion at a
constant engineering strain rate (again at $300\,\mathrm{K}$ and zero
pressure) to a target overall strain. We then return the simulation cell to an
approximately cubic shape via a biaxial strain in the other two directions.
Both stages of the compression last $0.5\,\mathrm{ns}$. By varying the target
strain (5%, 10%, 20%, 40%, 60% and 80%) we arrive at samples of material with
widely varying defect content with dislocation densities between zero and
$40\times 10^{16}\,\mathrm{m}^{-2}$ and stacking fault content from zero to
8%, calculated using the Ovito software 0965-0393-18-1-015012 . To further
vary the density and character of the defects we undertake brief annealing of
the deformed samples for $0.5\,\mathrm{ns}$ at a temperature of
$800\,\mathrm{K}$. Even this short anneal is long enough for the early stage
relaxation of the dislocation distributions to take place Race:2019ac .
Further details of our simulation approach are given in the companion paper
Race:ts .

Figure 1: Dislocation content of the cold-worked cells as a function of strain
for the as-worked (relaxed) and annealed ($0.5\,\mathrm{ns}$) cases.

The next step is to calculate a diffraction line profile from each of the
atomistic simulation cells. For an incoming x-ray beam with wavevector
$\mathbf{k}$ scattered into an outgoing beam with wavevector
$\mathbf{k}^{\prime}$ by an array of atoms with positions
$\\{\mathbf{r}_{i}\\}$ indexed by $i$, the relative phase difference for the
scattered wave from the $i^{\mathrm{th}}$ atom will be
$\exp(-\mathrm{i}\mathbf{Q}\cdot\mathbf{r}_{i})$, where
$\mathbf{Q}=\mathbf{k}^{\prime}-\mathbf{k}$ and $\hbar\mathbf{Q}$ is the
momentum transfer. Omitting a prefactor, which affects only the normalisation,
the amplitude of the scattered wave in the direction of $\mathbf{k}^{\prime}$
(assuming a fixed $\mathbf{k}$) is

$\Psi(\mathbf{Q})=\sum_{i}\mathrm{e}^{-\mathrm{i}\mathbf{Q}\cdot\mathbf{r}_{i}},$
(1)

assuming that all the atoms in our sample are of the same type and the angular
variation in scattered amplitude can be neglected. The scattered intensity is
given by
$I(\mathbf{Q})=\Psi^{\ast}(\mathbf{Q})\Psi(\mathbf{Q})=\sum_{i,j}\mathrm{e}^{-\mathrm{i}\mathbf{Q}\cdot\mathbf{r}_{ij}}$,
where $\mathbf{r}_{ij}=\mathbf{r}_{j}-\mathbf{r}_{i}$, and to arrive at the
equivalent of a powder pattern (or pattern for an untextured polycrystal) we
integrate this intensity over the full solid angle in $\mathbf{Q}$ assuming,
without loss of generality, that for each atom pair $\mathbf{r}_{ij}$ lies
along $\theta=0$.

This leads to the _Debye equation_ Debye:1915aa

$I(Q)=4\pi\sum_{i,j}\frac{\sin(Q\,r_{ij})}{Q\,r_{ij}},$ (2)

which allows us to calculate the whole line profile $I(Q)$ from a set of
atomic coordinates. Because $I(Q)$ depends on $\\{\mathbf{r}_{i}\\}$ only
through the interatomic separations $r_{ij}=|\mathbf{r}_{ij}|$, we can
implement (2) by calculating a histogram of the interatomic separations with
bin centres and frequencies $\\{(R_{s},f_{s})\\}$ so that the diffraction
profile is given by

$I(Q)=4\pi\sum_{s}f_{s}\,\frac{\sin(Q\,R_{s})}{Q\,R_{s}}\bigg{/}\sum_{s}f_{s}.$
(3)

Figure 2 shows part of the calculated line profiles for some of the simulated
materials in which the significant effect of the defect content on peak shape
is clearly visible.

Figure 2: Examples of line profiles for two low-index reflections in the
profiles calculated for the relaxed cold-worked material for a range of
strains.

We now have a set of simulation-based samples of deformed material for which
we know the exact defect content and have the corresponding diffraction line
profiles. All that remains is to apply line profile analysis and compare the
inferred defect content to the known, correct answer.

X-ray diffraction line profile analysis relies on the fact that the shapes of
the diffraction peaks for a sample of crystalline material depend on the size
of the crystallites and on the defects that they contain. The scattered
intensity for a given diffraction peak (indexed by $hkl$) is a convolution of
the effects of size and strain: $I_{hkl}=I_{hkl}^{\mathrm{size}}\ast
I_{hkl}^{\mathrm{strain}}$. These two effects can be separated because they
have a different dependence on diffraction order Warren:1952aa . The well
known Williamson-Hall method Williamson:1953aa  provides a simple analysis by
considering the full width at half maximum and integral breadth of each
diffraction peak, but thereby reduces the change in shape of the peaks to two
only parameters. Alternatively the Warren-Averbach method Warren:1950aa ;
Warren:1959aa  retains full information about peak shape by considering the
Fourier coefficients of the peaks
$A_{hkl}(L)=A_{hkl}^{\mathrm{size}}(L)\,A_{hkl}^{\mathrm{strain}}(L)$, where
$L$ denotes the Fourier variable. The strain broadening can be further
decomposed into contributions from different types of defect We can write, for
example,
$A_{hkl}^{\mathrm{strain}}(L)=A_{hkl}^{\mathrm{disl}}(L)\,A_{hkl}^{\mathrm{SF}}(L)$,
for contributions from dislocations and stacking faults. This decomposition is
possible because the Fourier coefficients for different defects vary in
different ways from reflection to reflection ($hkl$-dependence) and with
diffraction order.

A number of groups have implemented the Warren-Averbach method Delhez:1976wh ;
Berkum:1994ts ; Balzar:2004tg  and we make use of the approach named
_convolutional multiple whole profile_ (CMWP) analysis developed by Ungar and
co-workers Ungar:2001aa ; Ungar:1999aa ; Ungar:1984aa . CMWP uses physically
motivated profile functions for the Fourier coefficients $A_{hkl}(L)$ for
different defect types, parameterised in terms of their density and character
(e.g. burgers vector). The defect density is then inferred by optimising the
values of the parameters so that they give a theoretical line profile that
best matches the experimental profile. In the case of dislocations, Warren and
Averbach Warren:1950aa  showed that the effect of strain gives rise to Fourier
coefficients of the form

$A_{hkl}^{\mathrm{disl}}(L)\approx\exp(-G_{hkl}^{2}L^{2}\langle\varepsilon^{2}\rangle_{L}/2),$
(4)

where $G_{hkl}$ is the centre of the diffraction peak and
$\langle\varepsilon^{2}\rangle_{L}$ is the mean-square strain measured between
pairs of points in the material separated by a distance $L$ in a direction
perpendicular to the $(hkl)$ planes. This form for
$A_{hkl}^{\mathrm{disl}}(L)$ also assumes that we consider only small values
of $L$ and small strains. The derivation of (4) is rather involved and can be
found in Ref Warren:1950aa  in which (4) appears as Equation 14.

Wilkens Wilkens:1970aa  derived a form for the mean square strain for a so-
called restrictedly-random distribution of dislocations, in which dislocations
are randomly distributed within “cells” of material,

$\langle\varepsilon^{2}\rangle_{L}\approx\frac{\rho
Cb^{2}}{4\pi}\ln\left(\frac{R_{\mathrm{e}}}{L}\right),$ (5)

where $\rho$ is the dislocation density, $b$ is the length of the burgers
vector and $R_{\mathrm{e}}$ is the size of the “cell” for the restrictedly-
random distribution. $C$ is an average contrast factor, which quantifies the
strength of the effect of a given dislocation type on the diffraction peak
under study. The fundamental origins of this Krivoglaz-Wilkens strain function
in (5) are detailed in Refs. Groma:1997tk ; Groma:2016uc ; Zaiser:2001ur . We
investigate the interpretation of $R_{\mathrm{e}}$ further in our companion
paper Race:ts . The CMWP method is described in detail by Ungar et al.
Ungar:2001aa  and Ribarik et al. Ribarik:2020aa .

Figure 3 shows a comparison between the original diffraction profiles
generated from the simulation cells and the theoretical profiles corresponding
to the optimised values of the parameters in the model employed in the CMWP
process. In all cases, for a variety of strains, these fits are very good,
suggesting that there are no features peculiar to the profiles from the
simulated material that will prevent successful analysis via the CMWP process.

Figure 3: Comparison of original line profiles (thick, black lines) with CMWP
fits (thin, red lines) for samples annealed for $0.5\,\mathrm{ns}$ at several
strains. Profiles are offset vertically to aid readability. Note the log scale
in intensity.

We are now in a position to make a comparison between the defect content
inferred via the CMWP method and the true values, calculated from the
simulation cells. Figure 4(a) shows the comparison for the dislocation
densities for the relaxed and annealed material. Though quantitative agreement
is not perfect, the true dislocation density varies over two orders of
magnitude and CMWP does a good job of detecting this variation, generally
being within a factor of 2 of the true value.

(a)  (b)

Figure 4: A comparison of (a) dislocation densities and (b) stacking fault
fractions derived from line profile analysis with CMWP and the result of atom-
by-atom analysis for the relaxed (blue circles) and $0.5\,\mathrm{ns}$
annealed (orange squares) cells. Solid black lines join relaxed and annealed
points corresponding to the same initial strain to aid reading. The dashed
black lines indicate perfect agreement ($X_{\mathrm{Sim}}=X_{\mathrm{CMWP}}$)
between the atomistic analysis and CMWP and the dotted lines indicate a
difference of a factor of 2.

Figure 4(b) shows a comparison between the CMWP stacking fault fraction and
the fraction of atoms classified as having local hexagonal close packed
coordination in a common neighbour analysis. We see excellent agreement across
almost two orders of magnitude of stacking fault fraction.

By using atomistic simulations of cold-worked copper crystals to create
material for which we can exactly quantify the defect content, we are able to
directly test the performance of an important experimental method. We find
that the CMWP approach to line profile analysis (an implementation of the
Warren and Averbach approach) is able to simultaneously give an accurate
measure of both dislocation and stacking fault content across several orders
of magnitude of variation of defect density. Our results lend important new
support to the use of line profile analysis in the examination of deformed
crystalline material. Our companion paper Race:ts  contains further details of
the methods and analysis and investigates other aspects of the line profile
analysis, in particular probing the interpretation of the outer cut-off radius
($R_{\textrm{e}}$ in (5)) in the physical model for dislocation strain fields
used in the Warren and Averbach approach.

###### Acknowledgements.

C.P.R. was funded by a University Research Fellowship of The Royal Society.
T.U. was funded by the EPSRC Programme Grant (MIDAS: EP/S01702X/1). G.R. is
grateful for the support of OTKA grant K124926 funded by the Hungarian
National Research, Development and Innovation Office (NKFIH). Simulations were
carried out on the University of Manchester’s Computational Shared Facility.
The core research data, along with computer scripts implementing our methods,
are freely available for download Race:tb .

## References

  * [1] D. J. Bacon, D. M. Barnett, and R. O. Scattergood.  Anisotropic continuum theory of lattice defects.  Progress in Materials Science, 23:51–262, 1980.
  * [2] T. Ungar, H. Mughrabi, D. Rönnpagel, and M. Wilkens.  X-ray line-broadening study of the dislocation cell structure in deformed $[$001$]$-orientated copper single crystals.  Acta Metallurgica, 32(3):333–342, 1984.
  * [3] G. J. Ackland, G. Tichy, V. Vitek, and M. W. Finnis.  Simple n-body potentials for the noble metals and nickel.  Philosophical Magazine A, 56(6):735–756, 2013/01/25 1987.
  * [4] Steve Plimpton.  Fast parallel algorithms for short-range molecular dynamics.  Journal of Computational Physics, 117(1):1–19, 3 1995.
  * [5] Alexander Stukowski.  Visualization and analysis of atomistic simulation data with ovito–the open visualization tool.  Modelling and Simulation in Materials Science and Engineering, 18(1):015012, 2010.
  * [6] C P Race.  Atomistic simulations of grain boundary migration under recrystallisation conditions.  Modelling and Simulation in Materials Science and Engineering, 27(6):064002, 2019.
  * [7] C Race, T Ungar, and G Ribarik.  Validating x-ray line-profile defect analysis using atomistic models of deformed material.  Submitted to Physical Review Materials, 2022.
  * [8] P. Debye.  Zerstreuung von röntgenstrahlen.  Annalen der Physik, 351(6):809–823, 2021/02/04 1915.
  * [9] B. E. Warren and B. L. Averbach.  The separation of cold‐work distortion and particle size broadening in x‐ray patterns.  Journal of Applied Physics, 23(4):497–497, 2020/05/13 1952.
  * [10] G. K Williamson and W. H Hall.  X-ray line broadening from filed aluminium and wolfram.  Acta Metallurgica, 1(1):22–31, 1953.
  * [11] B. E. Warren and B. L. Averbach.  The effect of cold‐work distortion on x‐ray patterns.  Journal of Applied Physics, 21(6):595–599, 2020/05/14 1950.
  * [12] B. E. Warren.  X-ray studies of deformed metals.  Progress in Metal Physics, 8:147–202, 1959.
  * [13] R. Delhez and E. J. Mittemeijer.  The elimination of an approximation in the warren-averbach analysis.  Journal of Applied Crystallography, 9(3):233–234, 1976.
  * [14] J. G. M. van Berkum, A. C. Vermeulen, R. Delhez, T. H. de Keijser, and E. J. Mittemeijer.  Applicabilities of the warren-averbach analysis and an alternative analysis for separation of size and strain broadening.  Journal of Applied Crystallography, 27(3):345–357, 1994.
  * [15] D. Balzar, N. Audebrand, M. R. Daymond, A. Fitch, A. Hewat, J. I. Langford, A. Le Bail, D. Louer, O. Masson, C. N. McCowan, N. C. Popa, P. W. Stephens, and B. H. Toby.  Size-strain line-broadening analysis of the ceria round-robin sample.  Journal of Applied Crystallography, 37(6):911–924, 2004.
  * [16] T. Ungar, J. Gubicza, G. Ribarik, and A. Borbely.  Crystallite size distribution and dislocation structure determined by diffraction profile analysis: principles and practical application to cubic and hexagonal crystals.  J Appl. Cryst, 34(3):298–310, 2001.
  * [17] T. Ungar, I. Dragomir, A. Revesz, and A. Borbely.  The contrast factors of dislocations in cubic crystals: the dislocation model of strain anisotropy in practice.  Journal of Applied Crystallography, 32(5):992–1002, 1999.
  * [18] M. Wilkens.  The determination of density and distribution of dislocations in deformed single crystals from broadened x-ray diffraction profiles.  physica status solidi (a), 2(2):359–370, 2020/01/06 1970.
  * [19] I. Groma.  Link between the microscopic and mesoscopic length-scale description of the collective behavior of dislocations.  Physical Review B, 56(10):5807–5813, 09 1997.
  * [20] István Groma, Michael Zaiser, and Péter Dusán Ispánovity.  Dislocation patterning in a two-dimensional continuum theory of dislocations.  Physical Review B, 93(21):214110–, 06 2016.
  * [21] M. Zaiser, M. Carmen Miguel, and I. Groma.  Statistical dynamics of dislocation systems: The influence of dislocation-dislocation correlations.  Physical Review B, 64(22):224102–, 11 2001.
  * [22] Gábor Ribárik, Bertalan Jóni, and Tamás Ungár.  The convolutional multiple whole profile (cmwp) fitting method, a global optimization procedure for microstructure determination.  Crystals, 10(7), 2020.
  * [23] C P Race.  https://doi.org/10.5281/zenodo.6362428.