# TAX: Tendency-and-Assignment Explainer
for Semantic Segmentation with Multi-Annotators

Yuan-Chia Cheng
National Taiwan University
<EMAIL_ADDRESS>
Zu-Yun Shiau
National Taiwan University
<EMAIL_ADDRESS>
Fu-En Yang
National Taiwan University
<EMAIL_ADDRESS>
Yu-Chiang Frank Wang
National Taiwan University, NVIDIA
<EMAIL_ADDRESS>

###### Abstract

To understand how deep neural networks perform classification predictions,
recent research attention has been focusing on developing techniques to offer
desirable explanations. However, most existing methods cannot be easily
applied for semantic segmentation; moreover, they are not designed to offer
interpretability under the multi-annotator setting. Instead of viewing ground-
truth pixel-level labels annotated by a single annotator with consistent
labeling tendency, we aim at providing interpretable semantic segmentation and
answer two critical yet practical questions: “who” contributes to the
resulting segmentation, and “why” such an assignment is determined.

In this paper, we present a learning framework of Tendency-and-Assignment
Explainer (TAX), designed to offer interpretability at the annotator and
assignment levels. More specifically, we learn convolution kernel subsets for
modeling labeling tendencies of each type of annotation, while a prototype
bank is jointly observed to offer visual guidance for learning the above
kernels. For evaluation, we consider both synthetic and real-world datasets
with multi-annotators. We show that our TAX can be applied to state-of-the-art
network architectures with comparable performances, while segmentation
interpretability at both levels can be offered accordingly.

##  1 Introduction

Figure 1: Illustration of interpretable semantic segmentation with training
data annotated by multi-annotators. Practically, each annotator has his/her
labeling tendency. We aim to address semantic segmentation, while providing
interpretability at both annotator and assignment levels achieved (denoted in
blue and orange arrows, respectively).

Deep learning has presented its remarkable capability spanning a wide variety
of computer vision tasks, such as classification, segmentation, and object
detection. Despite the remarkable performances, deep learning models are often
viewed as unintelligible black boxes as the developers or the users are not
able to tell why and how such models function. In real-world applications, the
explainability of deep learning models and their prediction outputs are with
increasing demands. To this end, researchers begin to shift their efforts to
develop interpretable deep learning models.

Recent interpretable deep learning in computer vision has made notable
progress, particularly for classification models. A plethora of approaches
have been proposed to derive map-based [11, 30, 31, 38, 23, 14, 34, 4, 41, 20]
and example-based [6, 22, 21, 16, 5, 29] explanation, which points out how
local information is engaged in the image-level classification process. As for
interpretable segmentation, methods of [37, 33] can be viewed as map-based
manners built upon the backpropagation techniques or attention mechanisms.
However, most existing works towards interpretable learning models cannot
easily offer explanations to semantic segmentation with training image data
annotated by multiple and possibly diverse annotators.

In real-world applications, one typically requires multiple experts for data
labeling, since it would be time-consuming for one expert to annotate the
entirity of a dataset. Ideally, all the annotators are supposed to abide by a
standard annotating protocol or guideline (_e.g_., the PASCAL VOC labeling
guideline111http://host.robots.ox.ac.uk/pascal/VOC/voc2012/guidelines.html).
In practice, each annotator holds different labeling tendencies due to their
backgrounds. As noted in [39], some tend to over-segment image boundaries
while some might under-segment such data. Taking the multi-annotator context
into consideration, [39] proposes a coupled model comprised of an annotator
network estimating annotator confusion matrices when predicting the
segmentation output. However, existing methods like [39, 3] are not designed
to offer any type of interpretability during inference.

Instead of treating the label quality or labeling tendency equally the same
for each annotator, we aim at designing segmentation models with annotation
interpretability. To produce proper interpretations for multi-annotator
semantic segmentation as depicted in Figure 1, the segmentation model is
expected to answer two questions of interest: “who” contributes to the
resulting segmentation tendency (_i.e_., annotator-level explanation), and
“why” the model determines the corresponding tendency during inference
(_i.e_., assignment-level explanation). To realize the goal, we present a
novel learning framework of Tendency-and-Assignment Explainer (TAX), which
learns to perform interpretable segmentation, which jointly learns annotator-
specific masks and the corresponding convolution kernel subsets. More
specifically, the above kernel subsets are learned to describe the labeling
tendencies of each annotator. Moreover, with the annotator-dependent prototype
bank, our model allows to explain why particular annotator labeling tendency
is utilized during inference, realizing explainable and traceable segmentation
outputs. As confirmed by our experiments, our TAX can be applied to existing
segmentation models, producing satisfactory and explainable segmentation
outputs at annotator and assignment levels.

The contributions of this paper are highlighted below:

  * •

To the best of our knowledge, we are among the first to offer interpretability
for semantic segmentation models with the multi-annotator setting.

  * •

Instead of assuming ground-truth labels are with the same annotation style,
our TAX learning framework is designed to model specific labeling tendency,
aiming at providing both annotator and assignment-level explanations for
segmentation.

  * •

To offer interpretability at the annotator level, our TAX is able to learn
distinct convolutional kernel subsets given training data annotated by
different annotators.

  * •

Our TAX learns a prototype bank for describing image data, serving as
explainable guidance for the assignment and learning of the above kernel
subsets.

##  2 Related Works

###  2.1 Semantic Segmentation

A great number of deep learning approaches have been proposed for semantic
segmentation. As the pioneers, FCN [26] and U-Net [28] have demonstrated
excellent segmentation ability by adopting encoder-decoder architectures. To
better capture information at multiple resolutions/levels, PSPNet [40]
incorporates spatial pyramid pooling [25] at different grid scales, whereas
DeepLabv3 [9] utilizes Atrous Spatial Pyramid Pooling at different dilated
convolution rates. As the final extension of DeepLab models [7, 8, 9],
DeepLabv3+ [10] further exploits multi-scale contextual information by
designing a simplified decoder module which refines the segmentation results
along with object boundaries. On the other hand, EfficientPS [27] proposes a
semantic aggregation head with a Mask RCNN-based [18] head to encode and fuse
semantically rich multi-scale features.

In spite of impressive performances, existing models require collection of a
large number of ground-truth annotated data, which are often provided by
multiple users/experts. Expecting diverse segmentation tendencies, this leads
to the challenging multi-annotator semantic segmentation problem. Recently,
Probabilistic U-Net [24] and PHiSeg [3] exploit probabilistic CNN models to
model the inter-reader variations in segmentation labels. And [39] employs an
annotator network which estimates individual confusion matrices facilitating
the segmentation network to derive unobserved true label distribution.
Nevertheless, existing works generally focus on learning reliable information
across annotators, instead of providing outputs with proper explanation
(_e.g_., interpretable and traceable segmentation outputs). To the best of our
knowledge, no existing works can be directly applied to offer interpretability
under this multi-annotator semantic segmentation setting.

Figure 2: Overview of our Tendency-and-Assignment Explainer (TAX) learning
framework.  Given an input image $x_{i}$ annotated by $a_{i}$, the annotator
assigner $G$ predicts an annotator mask $M^{A}$ via the learned annotator-
dependent prototype bank $P$, guiding the segmentation model $F$ to learn
annotator-dependent convolution kernel set $C$ in its tendency learning layer
for providing explainable segmentation outputs. Best viewed in color.

###  2.2 Interpretable Deep Models

With the need of understanding the decision process of deep learning models, a
number of approaches have been developed. For image classification, the output
interpretability majorly falls into two types: map-based and example-based.
Map-based approaches [11, 30, 31, 38, 23, 14, 34, 4, 41, 20] offer
visualization maps to highlight the local parts of the input which contribute
to the classification result most. Specifically, some [30, 31, 38, 23, 14, 34,
4] estimate saliency maps by disturbing the input elements or backpropagating
the prediction to the input space, whereas some [41, 20] output attention maps
by incorporating CAM or attention mechanism into models. As for the example-
based branch [6, 22, 21, 16, 5, 29], one generally explains the prediction
with instance-wise examples in place of pixel-wise importance weights. For
instance, [5, 29] extract representative and discriminative part-wise
prototypes which are then taken as the references for conditional prediction.

Recent works like [37, 33] extend the above map-based techniques for
interpretable semantic segmentation. [37] applies guided backpropagation [32]
to FCN [26] for deriving attention maps, focusing on the object boundaries. On
the other hand, mitigating the time-inefficient post hoc overhead like
backpropagation, [33] explicitly fuses attentive modules into texture and
shape streams, where the texture stream handles coarse information for spatial
attention maps while the shape stream is in charge of fine boundary
information for shape attention maps. Nevertheless, the above methods are not
designed for multi-annotator semantic segmentation scenarios, and thus the
associated interpretability cannot be easily achieved.

##  3 Proposed Method

###  3.1 Problem Formulation and Method Overview

We first define the task of interpretable multi-annotator semantic
segmentation, including the notations used in this paper. In the multi-
annotator scenario, we have the $i$-th triplet as
$(x_{i},y_{i},a_{i})\in\\{X,Y,A\\}$ for training purposes. With a total of $N$
annotators available, each image $x_{i}$ is labeled by the corresponding
annotator $a_{i}=k$ ($k$ as the annotator index/label), providing the
associated ground-truth segmentation mask $y_{i}$. Note that the resolutions
of $x_{i}$ and $y_{i}$ are both of size $H\times W$ pixels.

In this paper, we propose a learning framework of Tendency-and-Assignment
Explainer (TAX), consisting of a tendency-preserving segmentation model $F$
and an annotator assigner $G$ (see Figure 2). We deploy and learn a set of
annotator-dependent convolution kernels $C$ in $F$, aiming to describe the
labeling tendency of each individual. On the other hand, annotator assigner
$G$ is to produce annotator masks $M^{A}$ which guide the learning of the
above kernels $C$. During inference, the annotator assigner $G$ predicts the
annotator mask $M^{A}$ for the input image, with $F$ inferring the
segmentation prediction using convolution kernels of the assigned annotators.
With the above inference process, interpretability at annotator and assignment
levels can be simultaneously achieved.

###  3.2 Learning to Describe Labeling Tendencies

To acquire tendency-explainable ability, we need the segmentation model to
capture individual labeling tendencies from multi-annotated training data. To
this end, we introduce a tendency-preserving segmentation model $F$ in Figure
2, which utilizes existing segmentation models (_e.g_., U-Net [28] or
DeepLabv3+ [10]) with only the last $3\times 3$ convolution layers modified.

In particular, we propose to learn the annotator-dependent convolution kernel
set $C$ in these final layers, since labeling tendencies are reflected at
finer feature resolutions. We have $C$ composed of annotator-dependent kernel
subsets $\\{C_{k}\\}_{k=1}^{N+1}$, with each kernel subset $C_{k}$ capturing
the $k$-th annotator’s labeling tendency. Note that the set size of $C$ is
$N+1$ rather than $N$, since we have the last kernel subset $C_{N+1}$
represent the labeling tendency shared by the annotators (_e.g_., background).

In the forwarding process of tendency learning layers, given the input feature
map $f$, we select different $C_{k}$ to perform convolution operation for each
pixel, and derive the output feature map $f^{\prime}\in\mathbbm{R}^{H\times
W\times D}$, where $D$ denotes the channel size. To select the proper kernel
subset for the forwarding process, our annotator assigner $G$ would produce an
annotator mask $M^{A}\in\mathbbm{R}^{H\times W}$, where each element
$m^{A}_{u,v}\in M^{A}$ indicates the annotator index for pixel $(u,v)$ (as
detailed in Sect. 3.3). If $m^{A}_{u,v}=k$ (recall that $k$ is the annotator
index/label), we then select the $k$-th kernel subset $C_{k}$ to perform
convolution operation for pixel $(u,v)$. Thus, at the tendency learning layer,
we have the output feature map calculated as:

$\displaystyle\mathbf{f^{\prime}_{u,v}}=conv(\mathbf{f_{u,v}},C_{k}\mid
m^{A}_{u,v}=k).$ (1)

To ensure that each kernel subset only captures the labeling tendency of the
corresponding annotator, we have $C_{k}$ only observe data labeled by the
exactly $k$-th annotator while $C_{N+1}$ learn from all the annotators. To
this end, we calculate the cross-entropy loss between the predicted
segmentation mask $\hat{y}_{i}$ and the ground-truth segmentation mask $y_{i}$
in an annotator-aware manner. More precisely, given training data
$(x_{i},y_{i},a_{i})$ (recall that one input image $x_{i}$ is annotated by one
single annotator $a_{i}$), we define the annotator-aware segmentation loss
$\mathcal{L}_{seg}^{(k)}$ for $k\in\\{1,2,\dotsb,N\\}$:

$\displaystyle\mathcal{L}_{seg}^{(k)}=\sum_{i}CE(\hat{y}_{i},y_{i})\mathbbm{1}_{k}(a_{i}=k).$
(2)

Note that $CE(\cdot,\cdot)$ denotes the cross-entropy loss, and
$\mathbbm{1}_{k}(\cdot)$ is an indicator function that outputs $1$ when the
image is annotated by the $k$-th annotator ($0$ otherwise). Thus, the
corresponding kernel subset $C_{k}$ is updated by:

$\displaystyle\theta_{C_{k}}\leftarrow\theta_{C_{k}}-\frac{\partial\mathcal{L}_{seg}^{(k)}}{\partial\theta_{C_{k}}},\quad
k\in\\{1,2,\dotsb,N\\}.$ (3)

On the other hand, since $C_{N+1}$ represents the kernel subset with tendency
shared by all annotators, we have the aggregated annotator-aware segmentation
loss $\mathcal{L}_{seg}$ as:

$\displaystyle\mathcal{L}_{seg}=\sum_{k=1}^{N}\mathcal{L}_{seg}^{(k)}\equiv\sum_{i}CE(\hat{y}_{i},y_{i}).$
(4)

Thus, this shared kernel subset $C_{N+1}$ is updated by:

$\displaystyle\theta_{C_{N+1}}\leftarrow\theta_{C_{N+1}}-\frac{\partial\mathcal{L}_{seg}}{\partial\theta_{C_{N+1}}}.$
(5)

As for updating the remaining part of the segmentation model $F$ (excluding
annotator-dependent kernels $C$), we apply the overall segmentation loss
(_i.e_., $\mathcal{L}_{seg}$) for optimizing the parameters
$\theta_{F\setminus C}$:

$\displaystyle\theta_{F\setminus C}\leftarrow\theta_{F\setminus
C}-\frac{\partial\mathcal{L}_{seg}}{\partial\theta_{F\setminus C}}.$ (6)

With the above learning strategies, each kernel subset $C_{k}\in C$ would
describe the desirable labeling tendency of each annotator, realizing
segmentation interpretability at the annotator level.

###  3.3 Learning to Assign Labeling Tendencies

As stated in Sect. 3.2, learning of annotator-dependent kernel subset $C_{k}$
is based on the assignment of the annotator mask $M^{A}$. We now describe how
we learn the annotator assigner $G$ for predicting this mask $M^{A}$, offering
segmentation interpretability at the assignment level.

##### Prototype-Bank Based Annotator Assigner $\boldsymbol{G}$

As depicted in Figure 2, we have an annotator assigner $G$ take the input
image $x_{i}$ and output the annotator mask $M^{A}\in\mathbbm{R}^{H\times W}$,
allowing the subsequent tendency learning layers to learn annotator-dependent
kernels. To offer assignment-traceable interpretability, this assigner $G$ is
based on the learning of annotator-dependent prototype bank $P$, which serves
as memory-bank-like [15] bases with visually explainable bases/prototypes.
Specifically, $P$ contains $N+1$ groups of annotator-dependent prototypes:

$\displaystyle
P=[P_{1},P_{2},\dotsb,P_{N+1}]\in\mathbbm{R}^{d\times((N+1)\times Q)},$ (7)

where prototype group $P_{k}\in P$ includes $Q$ trainable and explainable
prototypes $\mathbf{p}_{k}\in\mathbbm{R}^{d}$ describing local image patterns
associated with the tendency of the $k$-th annotator:

$\displaystyle
P_{k}=[\mathbf{p}_{k}^{1},\mathbf{p}_{k}^{2},\dotsb,\mathbf{p}_{k}^{Q}]\in\mathbbm{R}^{d\times
Q}.$ (8)

With the learned prototype bank $P$, we can derive $M^{A}$ for $x_{i}$,
reflecting the segmentation tendency for each local regions. To be more
precise, we feed $x_{i}$ into the prototype encoder $E$ to derive the
resulting feature map whose size is $h\times w\times d$, with each local
feature $\mathbf{E(x_{i})_{u,v}}\in\mathbbm{R}^{d}$ indicating the feature of
the location $(u,v)$. Then, we calculate the similarity between
$\mathbf{E(x_{i})_{u,v}}$ and the prototypes across annotators and identify
the annotator tendency of interest. That is, for each
$\mathbf{E(x_{i})_{u,v}}$, we calculate the cosine similarity score $s$ with
each prototype $\mathbf{p}_{k}\in P_{k}$, and select the highest score
$s_{k}^{\ast}$ for $k\in\\{1,2,\dotsb,N+1\\}$:

$\displaystyle s_{k}^{\ast}$ $\displaystyle=\mathop{\max}_{\mathbf{p}_{k}\in
P_{k}}sim(\mathbf{E(x_{i})_{u,v}},\>\mathbf{p}_{k}),$ (9)

where $sim(\cdot,\cdot)$ denotes the cosine similarity function. With $N+1$
selected scores, we thus have the soft version of the annotator mask
$M\in\mathbbm{R}^{h\times w\times(N+1)}$, in which each vector
$\mathbf{m}_{u,v}\in M$ indicates scores favoring the labeling tendencies
across $N+1$ annotators at the location $(u,v)$:

$\displaystyle\mathbf{m}_{u,v}=[s_{1}^{\ast},s_{2}^{\ast},\dotsb,s_{N+1}^{\ast}]\in\mathbbm{R}^{N+1}.$
(10)

To derive the final annotator mask $M^{A}\in\mathbbm{R}^{H\times W}$, we
simply apply the argmax operation on $M\in\mathbbm{R}^{h\times w\times(N+1)}$
for selecting the annotator index for each position and up-sample it to match
the image resolution size:

$\displaystyle M^{A}=UP_{r}(argmax(M)),$ (11)

where $argmax(\cdot)$ is performed along each vector $\mathbf{m}_{u,v}\in M$,
and $UP_{r}(\cdot)$ denotes the up-sampling operation with the ratio of
$r=(H/h,W/w)$.

To learn a reasonable $M^{A}$, we then introduce the pseudo annotator mask
$M^{pseu}\in\mathbbm{R}^{h\times w\times(N+1)}$ as the guidance, where
$M^{pseu}$ assigns a specific annotator for boundary prediction while
employing the shared annotator to handle the remaining consensus regions
(_e.g_., background). The motivation is that, the boundaries are the most
uncertain regions and thus are the major variance of segmentation tendencies
[39]. To produce $M^{pseu}$, we identify uncertain regions of an input image
$x_{i}$ by pre-training a vanilla segmentation model $F_{van}$ (_e.g_., a
plain U-Net [28]) on only $\\{X,Y\\}$ but $A$ (and then freeze for later
uses). Given an input $x_{i}$, the output of $F_{van}$ results in the
uncertainty mask $U\in\mathbbm{R}^{H\times W}$, where each element of $U$
indicates if the region is uncertain (represented by $1$) or not (represented
by $0$). After that, we assign specific labeling tendencies for those
uncertain regions while assigning the shared tendency for the remaining
consensus ones. Thus, we would assign the annotator index/label $a_{i}=k$ if
the element is $1$ on $U$, and specify $N+1$ otherwise, and then we convert
the annotator indexes to one-hot-encoded probabilities. Finally, we perform
down-sampling to obtain the final $M^{pseu}$ whose size is the same as that of
$M\in\mathbbm{R}^{h\times w\times(N+1)}$.

##### Updating Annotator Assigner $\boldsymbol{G}$

With the pseudo annotator mask $M^{pseu}$ as the guidance, we are able to
learn the annotator assigner $G$ using the following annotator assignment loss
$\mathcal{L}_{assign}$, which computes the cross-entropy loss ($CE$) between
$M_{i}$ and $M^{pseu}_{i}$ with respect to $x_{i}$:

$\displaystyle\mathcal{L}_{assign}=\sum_{i}CE(softmax(M_{i}),M^{pseu}_{i}).$
(12)

Note that $softmax(\cdot)$ transforms logits into a probability distribution
suitable for the cross-entropy loss calculation. Thus, the optimization of $G$
(only updating $E$ and $P$) can be formulated as follows:

$\displaystyle\theta_{E,P}\leftarrow\theta_{E,P}-\frac{\partial\mathcal{L}_{assign}}{\partial\theta_{E,P}}.$
(13)

It can be seen that, the motivation behind the above optimization is to
produce prototypes $\mathbf{p}_{k}\in P_{k}$ being substantially
representative to the uncertain local image regions of those labeled by the
$k$-th annotator ($k$ from $1$ to $N$), while prototypes $\mathbf{p}_{N+1}\in
P_{N+1}$ are encouraged to describe the remaining local regions (_e.g_.,
obvious foreground or background regions). With the completion of training
prototype banks $P$ and assigner $G$, our model would offer explainable and
traceable segmentation results at the assignment level.

###  3.4 Visual Interpretability during Inference

With learned $G$, $F$, and $C$, we now explain how interpretability at both
annotator and assignment levels can be achieved. For the annotator-level
explanation, the test input image is fed into $G$ and $F$ for deriving the
annotator mask and segmentation output. Thus, given a patch location of
interest, we are able to answer the question of “who” by pinning out the
annotator index on the annotator mask.

As for the assignment-level explanation, we are able to visualize the most
similar patches (via prototype identification with the patch size of
$\frac{H}{h}\times\frac{W}{w}$) for each pixel and answer the question of
“why”. That is, our model assigns the annotator for predicting the label for
that pixel, since the input patch is visually similar to those annotated by
the corresponding annotator during training.

##  4 Experiments

###  4.1 Datasets and Implementation Details

####  4.1.1 Datasets

PASCAL VOC 2012 [13] (referred to below as PASCAL VOC) contains 20 foreground
classes and one background class in realistic scenes. We adopt commonly used
augmented version [17] with 10,582/1,449/1,456 images for
training/validation/testing splits. We follow [10] and report the mean
intersection-over-union (mIoU) across the 21 classes on the validation split.

CVC-EndoSceneStill [36] (referred to below as EndoScene) is an endoluminal
scene segmentation benchmark for polyps from colonoscopy images, and contains
foreground (polyps) and background classes. The dataset is divided into
547/182/182 images for the training/validation/testing splits. We report the
mIoU across the two classes on the testing split.

LIDC-IDRI [1] holds 1018 lung computed tomography (CT) scans, resulting in
8,882/1,996/1,992 images as the training/validation/testing splits. Each image
is annotated by four thoracic radiologists, and we leave one label per image
for training data by following [39, 3], and we use single-label ground-truth
for evaluation [39]. We report the mIoU, mean-DICE (shortened as m-DICE) and
DICE (considering only the foreground class) on the testing split.

Since only the ground-truth labels of LIDC-IDRI are annotated by multiple
(four) experts, we explain how we manipulate multi-annotator labels for the
first two datasets. For PASCAL VOC and EndoScene, we partition the training
images into four subsets (each stands for one annotator), and we manipulate
their original masks with pre-defined labeling tendencies via four
morphological operations: dilated, eroded, simplified, and none, respectively.
As shown in Figure 3, dilated tends to over-segment, while eroded tends to
under-segment the boundaries. Simplified shows more straight contours, and
none indicates the original masks.

Figure 3: Examples of manipulated segmentation masks for imitating tendencies
of multi-annotators. Note that the foreground is shown in white for
simplicity.

####  4.1.2 Implementation Details

We implement our proposed tendency-preserving segmentation model $F$ with two
alternative architectures, U-Net [28] and DeepLabv3+ [10] (shortened as
DLv3+). For each architecture, we employ two encoder backbones pretrained on
ImageNet [12], ResNet101 [19] (shortened as Res101) and EfficientNet-B5 [35]
(shortened as Eff-B5). Prototype encoder $E$ and the vanilla segmentation
model $F_{van}$ are implemented simply with the U-Net [28] architecture with
ResNet50 [19] as the encoder backbone pretrained on ImageNet [12].

Figure 4: Case studies for interpretability on EndoScene and PASCAL VOC. Each
row presents an input image, segmentation prediction, and its interpretability
at annotator/assignment levels (_i.e_., explaining who/why). Note that the
original ground-truth masks are in red. Best viewed in color.

###  4.2 Case Studies for Interpretability

To verify the interpretability of our TAX, in Figure 4 we study annotator
(who) and assignment (why) explanations on EndoScene, PASCAL VOC, and LIDC-
IDRI.

##### EndoScene

As the first case shown in Figure 4(a), we see that our model produced a
simplified and rectangle-like segmentation mask, compared to the original,
non-manipulated ground-truth mask in red. One would raise the question about
whose labeling tendency results in such a segmentation output. After
examination, we notice that the predicted segmentation mask for this input
image was produced by the annotator with simplified labeling tendency (as
defined in Sect. 4.1.1). The annotator-specific results (_i.e_., annotator
masks) were as depicted in purple, reflecting his/her labeling tendency of
such annotation. Following such interpretability at the annotator level, the
subsequent question one would dig into is why the model chooses to utilize the
labeling tendency of this annotator. To answer this question, our annotator
assigner $G$ would be applied to identify the prototypes specifically selected
for performing such segmentation (see Sect. 3.3), suggesting the training
patches with high similarity to such input patches. This can be visually
verified by comparing the input and training patches in orange bounding boxes
(denoted with 1′, 2′), depicted in the first and last images/columns in Figure
4(a), respectively. Since EndoScene contains only polyps images, the selected
training patches do not exhibit significant visual appearance variety (_i.e_.,
limited assignment-level interpretability due to data). To better verify the
identifying ability of the annotator assigner, we then perform the case
studies on PASCAL VOC below.

##### PASCAL VOC

We now study a more challenging multi-class scenario on PASCAL VOC. In Figure
4(b), the second case study shows that our model predicted an over-segmented
result compared to the original, non-manipulated ground-truth segmentation
mask. The predicted annotator mask (in red) can be seen to be consistent with
the annotator with prior dilated labeling tendency (as defined in Sect.
4.1.1). Take the annotator mask around the helmet region for example. Our
annotator assigner identified the most similar training patch (a head with a
cap) via the learned annotator assigner and its prototype bank. Also, for the
motorcycle tire in the input image, its segmentation mask corresponds to the
same dilated labeling tendency, with the most similar training patch shown in
the last image of Figure 4(b).

Figure 5: Case studies for interpretability on LIDC-IDRI. Note that the
ground-truth masks are in red. Best viewed in color.

In the third case shown in Figure 4(c), the predicted segmentation exhibited
hybrid segmentation tendencies on the boundary, where the head of the horse
was segmented with fitting boundaries, while its feet were under-segmented.
For the annotator mask shown in Figure 4(c), we notice that these two parts
were indeed associated with different labeling tendencies, _i.e_., none and
eroded. As for tracing back the most similar training patches, the identified
patches are shown in the last two images of Figure 4(c), confirming the visual
and semantic similarities to the input ones. This case study confirms that,
our model can delicately assign different annotators for local segmentation
prediction. From the above case studies, we verify that our model not only
performs segmentation with interpretability at the annotator level, the
segmentation mask and its associated patch can be traced back to the training
ones with the same labeling tendency, realizing interpretability at the
assignment level.

##### LIDC-IDRI

We now perform case studies on the real-world multi-annotator dataset of LIDC-
IDRI. As shown in Figure 5, the predicted segmentation masks of the first and
second cases can be recognized as distinct tendencies, with annotator masks
shown in red and blue, respectively. Again, our model is able to explain such
segmentation tendency by tracing back to the most similar patches in the
training data. It is worth repeating that, prior works like [39, 3, 24] are
not designed to offer the above interpretability, and are designed to predict
a unified segmentation output without explanation of labeling tendency or
training data similarity. Thus, the use of our model would be desirable
especially for real-world applications such as medical image analysis.

###  4.3 Quantitative Analyses

####  4.3.1 Segmentation Performance

We now offer quantitative results to verify the effectiveness of our approach
using PASCAL VOC, EndoScene, and LIDC-IDRI. As noted in Sect. 3.2, our
approach is applicable to different segmentation architectures (_e.g_., U-Net
and DeepLabv3+), and thus we present segmentation results using different
backbones in Table 1. Particularly, baseline models in Table 1 denote non-
explainable vanilla segmentation models simply learned from the entire
training dataset with manipulated ground-truth labels. From Table 1, we see
that we achieved comparable mIoU as baseline models did, suggesting that
deploying our TAX learning strategy would not significantly affect the
segmentation performance, while interpretability can be additionally offered.

Furthermore, we conduct experiments on LIDC-IDRI for multi-annotator
segmentation. We compare with state-of-the-art methods of [24, 3, 2, 39],
which are designed for this setting yet not able to offer output
interpretability. As shown in Table 2, our model achieved comparable
segmentation performances as the competitors did. It is worth noting that, we
qualitatively and quantitatively verify the interpretability of our model in
Sections 4.2 and 4.3.2, while the above methods were not able to provide such
explanations. This again confirms the effectiveness of our method for real-
world semantic segmentation.

Archit. | Backbone | PASCAL VOC | EndoScene | LIDC-IDRI
---|---|---|---|---
Baseline | Ours | Baseline | Ours | Baseline | Ours
U-Net | Res101 | 69.28 | 69.31 | 82.80 | 82.60 | 54.02 | 54.52
Eff-B5 | 73.24 | 73.95 | 82.51 | 82.50 | 54.55 | 55.27
DLv3+ | Res101 | 74.86 | 74.72 | 81.90 | 82.80 | 55.20 | 54.93
Eff-B5 | 77.07 | 77.15 | 83.15 | 83.25 | 54.32 | 55.30

Table 1: Semantic segmentation performances on PASCAL VOC, EndoScene and LIDC-
IDRI in terms of mIoU (%). Note that our learning method can be applied to
different backbones and produce comparable performances, while additional
interpretability can be offered (as verified in Table 3).

Method | P. U-Net | PHiSeg | S. STAPLE | Zhang _et al_. | Ours
---|---|---|---|---|---
[24] | [3] | [2] | [39]
DICE | 52.38 | 54.08 | - | - | 53.25
m-DICE | - | - | 62.35 | 68.12 | 68.25

Table 2: Semantic segmentation on LIDC-IDRI in terms of DICE and m-DICE (%).
While our TAX is observed to achieve comparable results, the state-of-the-arts
for multi-annotator segmentation are not able to provide output explanation as
ours does.

####  4.3.2 Annotator and Assignment-Level Explanations

Earlier in Sect. 4.2, we qualitatively demonstrated that our method is capable
of offering visual explanations at the annotator and assignment levels. Here,
we assess the explanation abilities in quantitative manners.

To quantify the effectiveness of our annotator-level explanations, we examine
whether each individual convolution kernel subset $C_{k}$ exactly learns the
labeling tendency of the corresponding annotator $k$. In other words, we
choose to assess the segmentation performance of the learned tendency-
preserving segmentation model $F$ using only $C_{k}$. Specifically, for each
test image $x_{t}$ with ground truth $y_{t}$ annotated by a particular
labeling tendency (_i.e_., labeled by the annotator $k$), we then calculate
the mIoU between the predicted $\hat{y}_{t}$ and $y_{t}$. As presented in
Table 3, our method achieved superior performances by significant margins over
baseline models (_i.e_., vanilla segmentation models) which are not designed
to handle individual labeling tendency. This thus verifies that our tendency
learning layer succeeds in describing the labeling tendency of each annotator,
supporting our interpretability at the annotator level.

Archit. | Backbone | PASCAL VOC | EndoScene | LIDC-IDRI
---|---|---|---|---
Baseline | Ours | Baseline | Ours | Baseline | Ours
U-Net | Res101 | 62.60 | 66.01 | 77.84 | 80.00 | 53.20 | 53.42
Eff-B5 | 65.47 | 70.56 | 77.72 | 79.55 | 52.85 | 54.20
DLv3+ | Res101 | 67.23 | 69.15 | 77.53 | 80.00 | 53.85 | 53.70
Eff-B5 | 69.12 | 72.78 | 78.32 | 81.16 | 52.55 | 54.30

Table 3: Evaluation of segmentation with labeling tendency on PASCAL VOC,
EndoScene and LIDC-IDRI in mIoU (%). Note that, for PASCAL VOC and EndoScene,
pixel-level labels are manipulated by instructions presented in Sect. 4.1.1,
while LIDC-IDRI contains labels from multi-annotators.

As for the assessment of assignment-level explanations, we choose to evaluate
whether our learned annotator assigner $G$ is able to identify the visual
similarity between prototypes and input images. To conduct this quantitative
evaluation, we again have each test image $x_{t}$ with ground-truth labels
annotated by a particular labeling tendency (_e.g_., $a_{t}=k$, one of the
morphological operations). We then feed $x_{t}$ into $G$ to derive the
annotator mask $M^{A}$, followed by a majority vote across pixels $(u,v)$ in
$M^{A}$ for calculating the assignment accuracy. A considerably high
assignment accuracy of 85.80% was obtained on PASCAL VOC, confirming that our
annotator assigner is able to select proper prototypes visually similar to the
input patches. This thus further supports the interpretability of our model at
the assignment level.

##  5 Conclusion

In this paper, we aim at introducing interpretability to segmentation models
in the multi-annotator scenario, offering explanations to the segmentation
output at both annotator (who) and assignment (why) levels. With our proposed
framework of TAX, a proper set of convolution kernels would be derived for
describing the labeling tendency of each individual annotator, while the
learned prototype bank allows one to trace back the segmentation outputs to
the training image data. Our case studies on manipulated and real-world
datasets qualitatively supported the interpretability of our proposed method.
Moreover, our quantitative experiments verified the use of TAX with state-of-
the-art network models for producing satisfactory performances.

## References

  * [1] Samuel G Armato III, Geoffrey McLennan, Luc Bidaut, Michael F McNitt-Gray, Charles R Meyer, Anthony P Reeves, Binsheng Zhao, Denise R Aberle, Claudia I Henschke, Eric A Hoffman, et al.  The lung image database consortium (lidc) and image database resource initiative (idri): a completed reference database of lung nodules on ct scans.  Medical physics, 2011.
  * [2] Andrew J Asman and Bennett A Landman.  Formulating spatially varying performance in the statistical fusion framework.  IEEE transactions on medical imaging, 2012.
  * [3] Christian F Baumgartner, Kerem C Tezcan, Krishna Chaitanya, Andreas M Hötker, Urs J Muehlematter, Khoschy Schawkat, Anton S Becker, Olivio Donati, and Ender Konukoglu.  Phiseg: capturing uncertainty in medical image segmentation.  In MICCAI, 2019.
  * [4] Aditya Chattopadhay, Anirban Sarkar, Prantik Howlader, and Vineeth N Balasubramanian.  Grad-cam++: Generalized gradient-based visual explanations for deep convolutional networks.  In WACV, 2018.
  * [5] Chaofan Chen, Oscar Li, Chaofan Tao, Alina Jade Barnett, Jonathan Su, and Cynthia Rudin.  This looks like that: deep learning for interpretable image recognition.  In NeurIPS, 2019.
  * [6] Jianbo Chen, Le Song, Martin Wainwright, and Michael Jordan.  Learning to explain: An information-theoretic perspective on model interpretation.  In ICML, 2018.
  * [7] Liang-Chieh Chen, George Papandreou, Iasonas Kokkinos, Kevin Murphy, and Alan L Yuille.  Semantic image segmentation with deep convolutional nets and fully connected crfs.  In ICLR, 2015.
  * [8] Liang-Chieh Chen, George Papandreou, Iasonas Kokkinos, Kevin Murphy, and Alan L Yuille.  Deeplab: Semantic image segmentation with deep convolutional nets, atrous convolution, and fully connected crfs.  TPAMI, 2017.
  * [9] Liang-Chieh Chen, George Papandreou, Florian Schroff, and Hartwig Adam.  Rethinking atrous convolution for semantic image segmentation.  ArXiv, 2017.
  * [10] Liang-Chieh Chen, Yukun Zhu, George Papandreou, Florian Schroff, and Hartwig Adam.  Encoder-decoder with atrous separable convolution for semantic image segmentation.  In ECCV, 2018.
  * [11] Piotr Dabkowski and Yarin Gal.  Real time image saliency for black box classifiers.  In NeurIPS, 2017.
  * [12] Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei.  Imagenet: A large-scale hierarchical image database.  In CVPR, 2009.
  * [13] Mark Everingham, Luc Van Gool, Christopher KI Williams, John Winn, and Andrew Zisserman.  The pascal visual object classes (voc) challenge.  IJCV, 2010.
  * [14] Ruth C Fong and Andrea Vedaldi.  Interpretable explanations of black boxes by meaningful perturbation.  In ICCV, 2017.
  * [15] Dong Gong, Lingqiao Liu, Vuong Le, Budhaditya Saha, Moussa Reda Mansour, Svetha Venkatesh, and Anton van den Hengel.  Memorizing normality to detect anomaly: Memory-augmented deep autoencoder for unsupervised anomaly detection.  In ICCV, 2019.
  * [16] Sadaf Gulshad and Arnold Smeulders.  Explaining with counter visual attributes and examples.  In ICMR, 2020.
  * [17] Bharath Hariharan, Pablo Arbeláez, Lubomir Bourdev, Subhransu Maji, and Jitendra Malik.  Semantic contours from inverse detectors.  In ICCV, 2011.
  * [18] Kaiming He, Georgia Gkioxari, Piotr Dollár, and Ross Girshick.  Mask r-cnn.  In ICCV, 2017.
  * [19] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In CVPR, 2016.
  * [20] Zixuan Huang and Yin Li.  Interpretable and accurate fine-grained recognition via region grouping.  In CVPR, 2020.
  * [21] Jeya Vikranth Jeyakumar, Joseph Noor, Yu-Hsi Cheng, Luis Garcia, and Mani Srivastava.  How can i explain this to you? an empirical study of deep neural network explanation methods.  In NeurIPS, 2020.
  * [22] Atsushi Kanehira and Tatsuya Harada.  Learning to explain with complemental examples.  In CVPR, 2019.
  * [23] Andrei Kapishnikov, Tolga Bolukbasi, Fernanda Viégas, and Michael Terry.  Xrai: Better attributions through regions.  In ICCV, 2019.
  * [24] Simon AA Kohl, Bernardino Romera-Paredes, Clemens Meyer, Jeffrey De Fauw, Joseph R Ledsam, Klaus H Maier-Hein, SM Eslami, Danilo Jimenez Rezende, and Olaf Ronneberger.  A probabilistic u-net for segmentation of ambiguous images.  In NeurIPS, 2018.
  * [25] Svetlana Lazebnik, Cordelia Schmid, and Jean Ponce.  Beyond bags of features: spatial pyramid matching for recognizing natural scene categories.  In CVPR, 2006.
  * [26] Jonathan Long, Evan Shelhamer, and Trevor Darrell.  Fully convolutional networks for semantic segmentation.  In CVPR, 2015.
  * [27] Rohit Mohan and Abhinav Valada.  Efficientps: Efficient panoptic segmentation.  IJCV, 2021.
  * [28] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In MICCAI, 2015.
  * [29] Dawid Rymarczyk, Łukasz Struski, Jacek Tabor, and Bartosz Zieliński.  Protopshare: prototype sharing for interpretable image classification and similarity discovery.  ArXiv, 2020.
  * [30] Ramprasaath R Selvaraju, Michael Cogswell, Abhishek Das, Ramakrishna Vedantam, Devi Parikh, and Dhruv Batra.  Grad-cam: Visual explanations from deep networks via gradient-based localization.  In ICCV, 2017.
  * [31] Daniel Smilkov, Nikhil Thorat, Been Kim, Fernanda Viégas, and Martin Wattenberg.  Smoothgrad: removing noise by adding noise.  ArXiv, 2017.
  * [32] Jost Tobias Springenberg, Alexey Dosovitskiy, Thomas Brox, and Martin Riedmiller.  Striving for simplicity: the all convolutional net.  In ICLR workshop, 2015.
  * [33] Jesse Sun, Fatemeh Darbehani, Mark Zaidi, and Bo Wang.  Saunet: shape attentive u-net for interpretable medical image segmentation.  In MICCAI, 2020.
  * [34] Mukund Sundararajan, Ankur Taly, and Qiqi Yan.  Axiomatic attribution for deep networks.  In ICML, 2017.
  * [35] Mingxing Tan and Quoc Le.  Efficientnet: Rethinking model scaling for convolutional neural networks.  In ICML, 2019.
  * [36] David Vázquez, Jorge Bernal, F Javier Sánchez, Gloria Fernández-Esparrach, Antonio M López, Adriana Romero, Michal Drozdzal, and Aaron Courville.  A benchmark for endoluminal scene segmentation of colonoscopy images.  Journal of Healthcare Engineering, 2017.
  * [37] Kristoffer Wickstrøm, Michael Kampffmeyer, and Robert Jenssen.  Uncertainty modeling and interpretability in convolutional neural networks for polyp segmentation.  In MLSP, 2018.
  * [38] Shawn Xu, Subhashini Venugopalan, and Mukund Sundararajan.  Attribution in scale and space.  In CVPR, 2020.
  * [39] Le Zhang, Ryutaro Tanno, Mou-Cheng Xu, Chen Jin, Joseph Jacob, Olga Ciccarelli, Frederik Barkhof, and Daniel C Alexander.  Disentangling human error from the ground truth in segmentation of medical images.  In NeurIPS, 2020.
  * [40] Hengshuang Zhao, Jianping Shi, Xiaojuan Qi, Xiaogang Wang, and Jiaya Jia.  Pyramid scene parsing network.  In CVPR, 2017.
  * [41] Bolei Zhou, Aditya Khosla, Agata Lapedriza, Aude Oliva, and Antonio Torralba.  Learning deep features for discriminative localization.  In CVPR, 2016.


# Asymptotically Stable Non-Falling Solutions
of the Kapitza-Whitney Pendulum

Ivan Polekhin Electronic address<EMAIL_ADDRESS>Steklov Mathematical
Institute of Russian Academy of Sciences,
Moscow, Russia
Moscow Institute of Physics and Technology,
Dolgoprudny, Russia

## Abstract

The planar inverted pendulum with a vibrating pivot point in the presence of
an additional horizontal force field is studied. The horizontal force is not
assumed to be small or rapidly oscillating. We assume that the pivot point of
the pendulum rapidly oscillates in the vertical direction and the period of
these oscillations is commensurable with the period of horizontal force. This
system can be considered as a strong generalization of the Kapitza pendulum.
Previously it was shown that for any horizontal force there always exists a
non-falling periodic solution in the considered system. In particular, when
there is no horizontal force, this periodic solution is the vertical upward
position. In the paper we present analytical and numerical results concerning
the existence of asymptotically stable non-falling periodic solutions in the
system.

Keywords:  forced oscillation, the Kapitza pendulum, the Whitney pendulum,
stabilization, vibration.

##  1 Introduction

The problem of motion of a pendulum with a vibrating pivot point is one of the
well-studied non-linear dynamical systems. Like other classical systems — as
examples we can mention here the Duffing equation or the Van der Pol
oscillator — this system, firstly, is quite simple, which allows one to study
this system analytically, and secondly, some non-trivial dynamical effects can
be observed within the framework of this system.

Similarly to the the Duffing equation or the Van der Pol oscillator, the study
of the pendulum with a vibrating pivot point goes beyond the boundaries of
initial statements of the problems. Historically the problem of motion of a
pendulum with a vibrating base goes back to A. Stephenson [1], N.N. Bogolyubov
[2] and P.L. Kapitza [3, 4]. Let us also mention some recent papers on the
considering problem [5, 6, 7, 8]. The comprehensive bibliography can hardly be
listed here. However, one can find a relatively full overview of the papers on
the topic in [9, 10], including some references on the history of the problem.

In the presence of friction, for the pendulum with a vibrating base, one can
show that the vertical upward position (usually unstable) becomes
asymptotically stable, provided that the parameters of the system satisfy some
conditions. The proof follows from the classical method of averaging. If there
is no friction in the system, the proof is more complicated and also involves
KAM theory [11]. As a natural generalization of this system, one can consider
a planar mathematical pendulum with a vibrating point of suspension and in a
periodic in time horizontal force field.

If the horizontal force is rapidly oscillating, then we again can study this
system by means of the averaging theory. This problem has been considered
previously (see, for instance, [12]). In a more general statement of the
problem, it is not assumed that the horizontal force has a small period. In
[13] it was shown that for any horizontal periodic force, such that its period
is commensurate to the period of the vertical oscillation, there always exists
a periodic non-falling solution. When there is no horizontal force, then we
can consider the vertical upward equilibrium as a non-falling solution. Is it
always possible to make this non-falling solution stable by changing the
parameters of the vertical vibration? We will study this question in the
paper.

In the first section we present the equations of motion and some analytical
results on the problem. In particular, we will consider the case when the
horizontal force is weak. The case when the horizontal force is not weak is
considered numerically.

To be more precise, we consider the following two problems:

  1. 1.

Given any $2\pi$-periodic non-falling solution one can uniquely determine the
corresponding horizontal force, which guarantees the existence of this
solution. For which values of the amplitude of the vertical oscillation does
this solution become stable?

  2. 2.

Given any $2\pi$-periodic horizontal force. Again, we find numerically the
values of the amplitude such that there exists a periodic stable non-falling
solution.

In the conclusion section we will briefly discuss some conjectures related to
the problem.

##  2 Equations of motion

###  2.1 Planar pendulum

Let us consider a point of mass $m$ moving along a circle of radius $l$ in the
presence of a gravitational field. We assume that there is a force of viscous
friction acting on the point. If the circle is placed in a vertical plane,
then this system is the classical planar mathematical pendulum (Fig. 1). If
the pivot point of the pendulum moves along the vertical axis and the law of
motion is given by a function $h(t)$, then the position of the point is
defined as follows

$\displaystyle x=l\sin\varphi,$ $\displaystyle y=h(t)-l\cos\varphi.$

Here $\varphi$ is the angle between the vertical axis and the rod of the
pendulum. The kinetic and the potential energy ($T$ and $\Pi$, respectively)
of the system have the usual form

$\displaystyle
T=\frac{m}{2}(ml^{2}\dot{\varphi}^{2}+2\dot{h}\dot{\varphi}l\sin\varphi),\quad\Pi=-mgl\cos\varphi.$

Here $g$ is the acceleration of gravity. The usual Lagrangian equations of
motion for this system (with the Lagrangian function $L=T-\Pi$) have the
following form

$\displaystyle ml^{2}\ddot{\varphi}+m\ddot{h}l\sin\varphi=-mgl\sin\varphi.$

Figure 1: Inverted pendulum with a vibrating pivot point.

If we additionally consider viscous friction and a horizontal force acting on
the system, then the above equation have to be rewritten as follows

$\displaystyle
ml^{2}\ddot{\varphi}+m\ddot{h}(t)l\sin\varphi=-mgl\sin\varphi-\mu
l^{2}\dot{\varphi}+F(t)l\cos\varphi.$ (1)

Here $\mu$ is the coefficient of viscous friction and $F(t)$ is a function of
time. Without any loss of the generality, we can assume that $m=g=l=1$. Then
system (1) can be rewritten as follows

$\displaystyle\begin{split}&\dot{\varphi}=p-\dot{h}\sin\varphi,\\\
&\dot{p}=-\sin\varphi-\mu
p+\mu\dot{h}\sin\varphi+\dot{h}p\cos\varphi-\dot{h}^{2}\sin\varphi\cos\varphi+F\cos\varphi.\end{split}$
(2)

Below we assume that $h(t)$ is a rapidly oscillating function. To be more
precise, we assume that this function has the form

$h(t)=a\varepsilon\sin\frac{t}{\varepsilon}.$

For small $\varepsilon>0$ it is possible to consider the averaged system

$\displaystyle\begin{split}&\dot{\varphi}=p,\\\ &\dot{p}=-\sin\varphi-\mu
p-(a^{2}/4)\sin 2\varphi+F(s)\cos\varphi,\\\ &\dot{s}=1.\end{split}$ (3)

Here $s$ is an artificial time-like parameter (slow time) which is introduced
in order to distinguish two types of time-dependent functions.

##  3 Analytical results

###  3.1 Periodic solutions of the averaged system

In this section we obtain some analytical results on the existence of a
periodic solution (and its stability) for system (LABEL:eq2). First, we will
consider the same question for averaged system (LABEL:eq3). For instance, by
means of the classical Poincare method we will prove the existence of a
periodic asymptotically stable solution, provided that the perturbation is
small and non-autonomous. Therefore, from a result by N.N. Bogolyubov on the
existence of an integral manifold, we obtain the existence of a periodic
asymptotically stable solution in the original system.

We will also consider the case when the perturbation in (LABEL:eq3) is not
small. Here we will use a result by P. Torres on the existence and stability
of periodic solutions for the Duffing equation [14].

Everywhere below we assume that all functions are real analytic, except for
the cases when we explicitly assume that some function belongs to another
class of regularity.

######  Proposition 1.

Let function $F(t)$ in (LABEL:eq3) be as follows: $F(t)=cf(t)$, where $f$ is a
$2\pi$-periodic function. Let $\mu>0$ and $a^{2}>2$, then for sufficiently
small $c>0$, for system (LABEL:eq3), there exists a $2\pi$-periodic
asymptotically stable solution, and this solution coincides with the vertical
upward equilibrium for $c=0$.

###### Proof.

Let us denote the phase variable by $x=(\varphi,p)$, and the right hand side
by $G(t,x,c)$. The standard method of proving the existence of a periodic
solution for small $c$ is based on the application of the implicit function
theorem to the equation

$x(2\pi,x_{0};c)-x_{0}=0,$

here by $x(t,x_{0};c)$ we denote a the solution of the equation

$\dot{x}=G(t,x,c)$

with initial condition $x(0)=x_{0}$. For $c=0$ we have a trivial solution (the
vertical equilibrium of the pendulum). For small $c\neq 0$ there exists a
function $x_{0}(c)$ such that $x(2\pi,x_{0}(c),c)-x_{0}(c)=0$, provided that

$\det\left(\frac{\partial x(2\pi,x_{0};0)}{\partial x_{0}}-E\right)\neq 0.$

Here $E$ is the identity matrix. The determinant equals zero when the system
in variations (w.r.t. solution $\varphi=\pi$, $p=0$) has nontrivial periodic
solutions:

$\displaystyle\begin{split}&\dot{\varphi}=p,\\\
&\dot{p}=-\varphi\left(1+\frac{a^{2}}{2}\right)-\mu p.\end{split}$

It can be easily shown that for function

$H=\frac{p^{2}}{2}+\left(1+\frac{a^{2}}{2}\right)\frac{\varphi^{2}}{2}$

we almost everywhere have $\dot{H}=-\mu p^{2}<0$. Therefore, there exist no
non-trivial periodic solutions for the system in variations.

The obtained solutions will be asymptotically stable. Indeed, the vertical
position is asymptotically stable for $c=0$. If we consider this equilibrium
as a periodic solution, then we conclude that its characteristic multipliers
are strictly inside the unit circle. From the continuous dependence of the
characteristic multipliers on $c$, we obtain that for small $c$ the considered
solution is also asymptotically stable. ∎

When the friction is also small, i.e. if we put $c=\mu$ in (LABEL:eq3), then
it is also possible to prove the existence of an asymptotically stable
$2\pi$-periodic solution for small positive values of $\mu$.

This result follows from the following theorem proved by Lyapunov [15].

######  Theorem 1.

Let us have a non-autonomous Lyapunov system on a plane, i.e. a system of
differential equations of the following form

$\displaystyle\begin{split}&\dot{x}=-\lambda y+f_{x}(x,y)+\mu
F_{x}(t,x,y,\mu),\\\ &\dot{y}=\lambda x+f_{y}(x,y)+\mu
F_{y}(t,x,y,\mu).\end{split}$ (4)

Here $F_{x},F_{y}$ are $2\pi$-periodic functions in $t$. Functions $f_{x}$ и
$f_{y}$ are of order no less than $two$ in $x,y$. We also assume that

$f_{x}=-\frac{\partial S}{\partial y},\quad f_{y}=\frac{\partial S}{\partial
x}$

If $\lambda\not\in\mathbb{Z}$, then there exists a unique $2\pi$-periodic
solution of system (LABEL:eq7). This solution depends analytically on $\mu$
for small $\mu$, and for $\mu=0$ we obtain the trivial solution $x=y=0$.
Moreover, if

$\mu\int\limits_{0}^{2\pi}\left.\left(\frac{\partial f}{\partial
x}+\frac{\partial F}{\partial y}\right)\right|_{x=y=\mu=0}\,dt<0,$

then all solutions in the obtained one-parameter family of periodic solutions
are asymptotically stable.

If we apply this theorem to system (LABEL:eq3), we obtain the following
proposition.

######  Proposition 2.

Let $(1+a^{2}/2)^{1/2}\notin\mathbb{Z}$, $F(t)=\mu f(t)$, where $f(t)$ is a
$2\pi$-periodic function. Then for small $\mu>0$ there is a $2\pi$-periodic
asymptotically stable solution of (LABEL:eq3).

###### Proof.

System (LABEL:eq3) can be presented in the following form

$\displaystyle\begin{split}&\dot{x}=-(1+a^{2}/2)^{1/2}y+X(y)-\mu y+\mu
f(t)\cos(y/(1+a^{2}/2)^{1/2}),\\\ &\dot{y}=(1+a^{2}/2)^{1/2}x.\end{split}$

Where $X(y)$ is an analytic function of order no less than two in $y$. If
$\mu=0$, then this system is Hamiltonian. From Theorem 1 we obtain the
existence of a periodic solution. Since $\mu>0$, then

$-2\pi\mu<0.$

Therefore, this solution is asymptotically stable. ∎

Let us now consider the case when it is not assumed that $F(t)$ is small. We
will use the following version of a result by P. Torres [14] (see also [16]).

Given a Duffing equation

$\displaystyle\ddot{x}+c\dot{x}+g(t,x)=0,$ (5)

where $c>0$, $g\colon\mathbb{R}^{2}\to\mathbb{R}$ is $2\pi$-periodic in $t$,
by $\Omega_{k,c}$ we denote the following set

$\Omega_{k,c}=\left\\{f\in L^{k}(0,2\pi)\colon f\succ
0,\|f\|_{k}<\left(1+\frac{c^{2}}{4}\right)K\left(\frac{2k}{k-1}\right)\right\\}.$

Here $f\succ 0$ means that $f\geqslant 0$ (a.e.) and $f>0$ in a subset of
positive measure.

$K(q)=\begin{cases}\frac{1}{q\cdot(2\pi)^{2/q}}\left(\frac{2}{2+q}\right)^{1-2/q}\left(\frac{\Gamma(1/2)}{\Gamma(1/2+1/q)}\right)^{2},&\mbox{
if }1\leqslant q<+\infty,\\\ 2/\pi,&\mbox{ if }q=+\infty.\end{cases}$

######  Theorem 2.

Let $\alpha>\beta$ be two numbers such that $g(t,\alpha)<0$ and $g(t,\beta)>0$
for all $t$ and for some $1\leqslant k\leqslant+\infty$ there exists
$f\in\Omega_{k,c}$ and for all $x\in[\beta,\alpha]$

$\frac{\partial g(t,x)}{\partial x}\leqslant f(t)\mbox{ {(a.e.).} }$

Then system (5) has at least one asymptotically stable solution $x(t)$ and
$x(t)\in(\beta,\alpha)$ for all $t$.

As a corollary from this theorem we obtain the following result on the
existence of an asymptotically stable $2\pi$-periodic solution of (LABEL:eq3).
First, let us introduce the following notation

$\Phi(\varphi)=-\sin\varphi-\frac{a^{2}}{4}\sin 2\varphi.$

If $a^{2}>2$, then function $\Phi$ has two local maxima inside the interval
$(\pi/2,3\pi/2)$ and $\Phi>0$ at these points; $\Phi$ also has two local
minima where $f<0$. By $\lambda_{1}$ and $\lambda_{2}$ we denote the following
numbers

$\lambda_{1}=\frac{-1+\sqrt{1+2a^{4}}}{2a^{2}},\quad\lambda_{2}=\frac{-1-\sqrt{1+2a^{4}}}{2a^{2}}.$

For $a^{2}\in(0,2)$ two critical points (inside $[0,2\pi]$) of $\Phi(\varphi)$
are as follows

$\varphi_{min}^{1}=\mathrm{arccos}(\lambda_{1}),\quad\varphi_{max}^{1}=2\pi-\mathrm{arccos}(\lambda_{1}).$

As $a^{2}$ tends to $0$, value $\varphi_{min}^{1}$ tends to $\pi/2$ and
$\varphi_{max}^{1}$ tends to $3\pi/2$. As $a^{2}$ tends to $\infty$,
$\varphi_{min}^{1}$ tends to $\pi/4$ and $\varphi_{max}^{1}$ tends to
$7\pi/4$. If $a^{2}>2$, then we have two more additional critical points

$\varphi_{max}^{2}=\mathrm{arccos}(\lambda_{2}),\quad\varphi_{min}^{2}=2\pi-\mathrm{arccos}(\lambda_{2}).$

As $a^{2}$ tends to $2$, $\varphi_{min}^{2}$ tends to $\pi$ and
$\varphi_{max}^{2}$ also tends to $\pi$. As $a^{2}$ tends to $\infty$,
$\varphi_{min}^{2}$ tends to $5\pi/4$ and $\varphi_{max}^{2}$ tends $3\pi/4$.

######  Proposition 3.

For any $2\pi$-periodic $F(t)$ such that $|F(t)|\leqslant 2/\pi$ and
$-F(t)\cos\varphi^{2}_{max}<\Phi(\varphi^{2}_{max})$,
$\Phi(\varphi^{2}_{min})<-F(t)\cos\varphi^{2}_{min}$ for all $t$, there exists
a $2\pi$-periodic non-falling asymptotically stable solution of (LABEL:eq3).

###### Proof.

Let us put $\alpha=\varphi^{2}_{min}$ and $\beta=\varphi^{2}_{max}$. Then
inequalities $g(t,\alpha)<0$ and $g(t,\beta)>0$ hold. We have

$g(t,\varphi)=\sin\varphi+\frac{a^{2}}{4}\sin 2\varphi-F(t)\cos\varphi,$

and

$\frac{\partial g}{\partial\varphi}=\cos\varphi+\frac{a^{2}}{2}\cos
2\varphi+F(t)\sin\varphi.$

We also have that

$\cos\varphi+\frac{a^{2}}{2}\cos 2\varphi\leqslant 0$

for all $\varphi\in[\varphi^{2}_{max},\varphi^{2}_{min}]$. Finally, we can put
$f(t)=2/\pi$, $k=+\infty$ and apply Theorem 2. ∎

###  3.2 Periodic solutions of the original system

Let us now discuss the connection between the existence of periodic solutions
of the averaged (LABEL:eq3) and the original (LABEL:eq2) systems. The basic
result which establish this connection is a classical result by N.N.
Bogolyubov on the existence of an integral manifold. A more detailed
explanation of the results presented below can be found in [17, 18].

Given a system of ordinary differential equations in the standard form of the
method of averaging

$\frac{dx}{dt}=\varepsilon X(t,x),$ (6)

where $x\in\mathbb{R}^{n}$, $t\in\mathbb{R}$, $\varepsilon$ is a parameter,
$X\colon\mathbb{R}\times\mathbb{R}^{n}\to\mathbb{R}^{n}$. Let function $X$ be
a $2\pi$-periodic in $t$. We call system (6) the original system. Let us also
consider the averaged system

$\frac{dx}{d\tau}=X_{0}(x),$ (7)

where function $X_{0}$ is the average of function $X$ w.r.t. the periodic
variable $t$ and $\tau=\varepsilon t$

$X_{0}(x)=\frac{1}{2\pi}\int\limits_{0}^{2\pi}X(t,x)\,dt.$

Let us assume that averaged system (7) has a periodic solution
$\tilde{x}(\tau)$ of period $2\pi$ and the characteristic multipliers of this
solution are strictly inside the unit circle. Then for sufficiently small non-
negative $\varepsilon$ there exists an integral manifold close to solution
$\tilde{x}$. Let us recall that we say that a manifold is integral if any
solution starting at this manifold never leaves it. In coordinates, this
integral manifold, which we denote by $S$, can be presented as follows

$x=f_{\varepsilon}(t,\theta),$

where $f_{\varepsilon}$ is $2\pi$-periodic in both variables. The period in
$t$ is equal to the corresponding period of function $X(t,x)$, $\theta$ is an
angular variable and we can assume without loss of the generality that
function $f_{\varepsilon}(t,\theta)$ is $2\pi$-periodic in $\theta$.

In a sufficiently small neighborhood of the periodic solution our integral
manifold is unique. Moreover, this integral manifold is asymptotically stable,
i.e. any solution of the original system asymptotically tends to $S$, provided
that the initial conditions for this solution are located sufficiently close
to $S$. At the same time, in the general case, we cannot conclude that all
solutions on $S$ are periodic for arbitrary $\varepsilon$.

However, for system (LABEL:eq2) it can be proved that integral manifold $S$ is
stratified into $2\pi$-periodic solutions, provided that $\varepsilon=1/k$ and
$k\in\mathbb{N}$ is a sufficiently large number. Note, that in our case
angular variable $\theta$ corresponds to variable $s$, which defines the
initial moment of time for our original system. Since the original and the
averaged equations for $s$ coincides, then from the asymptotic stability of
$S$ we obtain the asymptotic stability for every periodic solution on $S$
(provided that we do not change the initial condition for variable $s$).

In other words, if the periods of the vertical oscillation and the horizontal
force are commensurable and the averaged system has an asymptotically stable
$2\pi$-periodic solution, then the original system also has an asymptotically
stable $2\pi$-periodic solution, provided that $k$ is large.

##  4 Numerical results

###  4.1 Stability of a given periodic motion

The above analytical results hold only in the case when the external
horizontal force is weak or some additional assumptions hold. Below we present
numerical results for the cases when neither $\mu$ nor $F(t)$ in system
(LABEL:eq3) is small.

By a direct calculation, we obtain the following result showing that any given
periodic non-falling solution can be obtained as a solution of our system for
some periodic horizontal time-dependent force.

######  Lemma 1.

For any $2\pi$-periodic function $\varphi(t)$ such that
$\varphi(t)\in(\frac{\pi}{2},\frac{3\pi}{2})$ for all $t$, there exists a
unique $2\pi$-periodic function $F(t)$ such that $\varphi(t)$ is a solution of
(LABEL:eq3). Function $F(t)$ has the following form:

$F(t)=\frac{\ddot{\varphi}(t)+\sin\varphi(t)+\mu\dot{\varphi}(t)+(a^{2}/4)\sin
2\varphi(t)}{\cos\varphi(t)}.$

Let us assume that the values of parameters $\mu$ and $a$ in system
(LABEL:eq3) are given. Given any non-falling $2\pi$-periodic trajectory of
motion $\varphi(t)$ ($\varphi(t)\in(\pi/2,3\pi/2)$ for all $t$), we find
$F(t)$ such that $\varphi(t)$ is a solution of (LABEL:eq3). Here we put
$p(t)=\dot{\varphi}(t)$. For instance, if we put
$\varphi(t)=\pi-\pi/4\cdot\cos t$, then we consider a periodic solution of
amplitude $\pi/4$. For some values of $a$ our solution may become
asymptotically stable. In this case, in the original system we also have a
non-falling $2\pi$-periodic solution and these two solutions of systems
(LABEL:eq2) and (LABEL:eq3) are close provided that $\varepsilon=1/k$ is
small.

(a) $\mu=5$

(b) $\mu=3$

(c) $\mu=1$

(d) $\mu=0.1$

Figure 2: Regions of asymptotic stability for various $\mu$.

We will consider the following one-parameter family of solutions:

$\varphi_{1}(t)=\pi-A\cos t,\quad A\in[0,\pi/2).$

At the same time, we will consider various values of parameter $\mu>0$. In
Figure 2 one can find the regions of asymptotic stability for this problem.

Note that some points of the region of asymptotic stability are arbitrarily
close to vertical line $A=\pi/2$. It means that for any amplitude $A<\pi/2$ we
can stabilize the corresponding solution by choosing some $a$ from a bounded
interval.

Moreover, the distance between the consecutive points in which the region of
asymptotic stability is arbitrarily close to line $A=\pi/2$ equals $2\pi$.
Based on the numerical results, it is also possible to conclude that for the
first point in this sequence we have $a\to\pi$ as $\mu\to 0$.

###  4.2 Stability of a periodic non-falling solution for a given horizontal
force

The second problem that we study numerically can be formulated as follows. For
any $\mu>0$ и $F\equiv 0$ there exists the vertical equilibrium of system
(LABEL:eq3) for which $\varphi=\pi$. This equilibrium is asymptotically stable
when $a>\sqrt{2}$. It can be shown that for any $T$-periodic horizontal force
$F$ there exists at least one non-falling $T$-periodic solution. If $F\equiv
0$ then the vertical equilibrium can be considered as the required periodic
non-falling solution solution. We will numerically find the values of $a$ such
that the corresponding periodic solution is asymptotically stable: we will
start from know condition $a>\sqrt{2}$ (for $A=0$) and will numerically
continue this condition for larger $A$.

To be more precise, if our force $F$ depends on some parameter (denoted by
$A$), then we will obtain a plot of $a$ as a function of $A$. For any given
$A$, the corresponding value of $a$ can be described as follows: if the
amplitude of the vertical oscillations is slightly less than $a$, then the
corresponding periodic solution is non-stable, if the amplitude of the
vertical oscillations is slightly greater than $a$, then the corresponding
periodic solution is asymptotically stable. Everywhere below we assume that
$T=2\pi$.

Let us present a more detailed explanation of the algorithm. For any given
force $F$ we have to find a periodic non-falling solution of system
(LABEL:eq3). Let us consider the following function
$\Phi\colon\mathbb{R}\times\mathbb{R}\to\mathbb{R}$

$\Phi(\varphi_{0},p_{0})=\left((\varphi(2\pi;\varphi_{0},p_{0})-\varphi_{0})^{2}+(p(2\pi;\varphi_{0},p_{0})-p_{0})^{2}\right)^{1/2}.$

Here $\varphi(t;\varphi_{0},p_{0})$ are $p(t;\varphi_{0},p_{0})$ the
components of the solution of system (LABEL:eq3) when $\varphi(0)=\varphi_{0}$
и $p(0)=p_{0}$. The periodic solutions correspond to the zeros of $\Phi$. We
will use the method of numerical continuation and will try to find zeros close
to the zeros obtained for close values of $A$. In other words, we restrict the
search for periodic solutions to a neighborhood of the periodic solution
obtained on the previous step (for the previous value of $A$).

When the zeros are found, we obtain numerically the monodromy matrices of the
corresponding periodic solutions and check whether or not all characteristic
multipliers belong to the unit circle. The required values of $a$ is
determined by means of the bisection method.

Two level set plots of two different functions $\Phi$ and their zeros are
shown in Figure 3.

(a) $A$ is small, $a$ is large ($\mu=1$, $A=0.1$, $a=2$).

(b) $a$ is small, $A$ is large ($\mu=1$, $A=2$, $a=1$).

Figure 3: Asymptotically stable (yellow) and unstable (blue) $2\pi$-periodic
solutions and the corresponding zeros of function $\Phi$.

The resulting plots are shown in Figure 4. In Figure 5 one can find how an
unstable periodic solution bifurcates and two unstable and one asymptotically
stable periodic solutions appear as we cross the curve presented in Figure 4
for $\mu=0.1$.

The problem of finding of all periodic non-falling solutions for a given $A$
is more difficult from the computational point of view. Therefore, we mainly
focus on the numerical continuation of a known solution. The following lemma
allows one to simplify the search of a periodic solution for system
$\eqref{eq3}$. This Lemma has been used to obtain Figure 5.

Figure 4: $a$ vs. $A$ plots for $\mu$ equals $0.1$, $1$, and $3$. The regions
above the plots corresponds to the values of $a$ and $A$ such that there
exists at least one $2\pi$-periodic asymptotically stable non-falling
solution.

######  Lemma 2.

A periodic non-falling solution of system (LABEL:eq3) can exist only inside
the rectangle $\varphi\in[\pi/2,3\pi/2]$, $p\in[-P,P]$, where

$P=\frac{1+a^{2}/4+\max|F|}{\mu}.$

###### Proof.

The first condition is the definition of a non-falling solution. Let us show
that a periodic solution cannot exist in the region where $p>P$. Indeed, in
all points of this region we have $\dot{p}<0$ (it follows from system
(LABEL:eq3)). Therefore, we obtain that the solution cannot return to its
initial condition. Similarly one can consider the case when $p<-P$. ∎

(a) $a=1.4220$

(b) $a=1.4240$

Figure 5: Bifurcation of periodic solutions for $A=0.1$ and $\mu=0.1$.

##  5 Conclusion

A generalization of the well known Kapitza pendulum has been considered. In
the case when there is no additional horizontal force the considered and
system is obviously equivalent to the Kapitza pendulum. Some results on the
existence of asymptotically stable periodic non-falling solutions are
presented. In particular, the obtained numerical results (Figure 4) generalize
classical conditions for the asymptotical stability of the upward vertical
equilibrium of the Kapitza pendulum.

However, the considered system is far from being completely studied. In the
conclusion section we would like to outline some open problem, which will be
studied elsewhere.

First, it would be interesting to prove that for any periodic (and, therefore,
bounded) force $F(t)$ there exists a sufficiently large $a$ such that system
(LABEL:eq3) has a periodic non-falling solution provided that
$\varepsilon=1/k$ is small. It is interesting to find the correspondence
between the norm of $F(t)$ and the minimal possible value of $a$ such that
this solution exists.

Second, it may be of some practical interest to prove results concerning the
qualitative view of plots shown in Figure 4. Is it true that these plots are
always monotonously increasing? Is it possible to find several separated
regions of stability for a given $A$ (similarly to the cases shown in Figure
2)? It will be interesting to compare Figure 2 with the well known Ince-Strutt
diagram for the Mathieu equation (see, for instance, [19]).

Third, it is known that under some additional assumptions, there exist at
least two periodic solutions of the Kapitza-Whitney pendulum [13]. Is it true
that the existence of two unstable periodic non-falling solutions guarantees
the existence of at least one stable periodic non-falling solution? Note that
the mentioned two solutions are unstable if there is no friction in the system
[20]. From the above results it follows that this statement holds if $F(t)$ is
small.

As a last open problem in our short list we would like to mention the problem
of existence of so-called subharmonic solutions (the solution is called
subharmonic if this solution is $nT$-periodic and not $kT$-periodic for all
$k<n$, $k,n\in\mathbb{N}$). It would be interesting to compare the minimal
values of $a$ which stabilize the solution for $T$-periodic and $nT$-periodic
solutions.

Finally, we can add that many questions that can be asked about the Kapitza
pendulum can be easily reformulated for the Kapitza-Whitney pendulum.

## Acknowledgements

This work was supported by the Russian Science Foundation (Project no.
19-71-30012).

## References

  * [1] Andrew Stephenson  ‘‘On induced stability’’  In _The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science_ 15.86 Taylor & Francis, 1908, pp. 233–236
  * [2] N.. Bogolyubov  ‘‘Perturbation theory in nonlinear mechanics’’  In _Collection of Papers of Inst. Construction Mekhaniki Akad. Nauk UkrSSR_ 14, 1950, pp. 9–34
  * [3] P.. Kapitsa  ‘‘Dynamic stability of the pendulum when the point of suspension is oscillating’’  In _Sov. Phys. JETP_ 21, 1951, pp. 588
  * [4] P.. Kapitsa  ‘‘The pendulum with an oscillating pivot point’’  In _Uspekhi fizicheskikh nauk_ 44.7, 1951, pp. 7–20
  * [5] Gerson Cruz Araujo and Hildeberto E Cabral  ‘‘Parametric Stability of a Charged Pendulum with an Oscillating Suspension Point’’  In _Regular and Chaotic Dynamics_ 26.1 Springer, 2021, pp. 39–60
  * [6] Zvi Artstein  ‘‘The pendulum under vibrations revisited’’  In _Nonlinearity_ 34.1 IOP Publishing, 2021, pp. 394
  * [7] A.. Belyaev et al.  ‘‘Classical Kapitsa’s problem of stability of an inverted pendulum and some generalizations’’  In _Acta Mechanica_ Springer, 2021, pp. 1–17
  * [8] Hildeberto E Cabral and Adecarlos C Carvalho  ‘‘Parametric stability of a charged pendulum with oscillating suspension point’’  In _Journal of Differential Equations_ 284 Elsevier, 2021, pp. 23–38
  * [9] Eugene I Butikov  ‘‘On the dynamic stabilization of an inverted pendulum’’  In _American Journal of Physics_ 69.7 American Association of Physics Teachers, 2001, pp. 755–768
  * [10] Anatoly Mykhailovych Samoilenko  ‘‘N.N. Bogolyubov and non-linear mechanics’’  In _Russian Mathematical Surveys_ 49.5 IOP Publishing, 1994, pp. 109
  * [11] BS Bardin and AP Markeyev  ‘‘The stability of the equilibrium of a pendulum for vertical oscillations of the point of suspension’’  In _Journal of Applied Mathematics and Mechanics_ 59.6 Elsevier, 1995, pp. 879–886
  * [12] Vladimir Burd  ‘‘Method of averaging for differential equations on an infinite interval: theory and applications’’  ChapmanHall/CRC, 2007
  * [13] Ivan Yu Polekhin  ‘‘The Method of Averaging for the Kapitza–Whitney Pendulum’’  In _Regular and Chaotic Dynamics_ 25.4 Springer, 2020, pp. 401–410
  * [14] Pedro J Torres  ‘‘Existence and stability of periodic solutions of a Duffing equation by using a new maximum principle’’  In _Mediterranean Journal of Mathematics_ 1.4 Springer, 2004, pp. 479–486
  * [15] I G Malkin  ‘‘Methods of Lyapunov and Poincaré in the Theory of Nonlinear Vibrations (in Russian)’’  Moscow–Leningrad: Gostekhizdat, 1949
  * [16] Franic Ikechukwu Njoku and Pierpaolo Omari  ‘‘Stability properties of periodic solutions of a Duffing equation in the presence of lower and upper solutions’’  In _Applied mathematics and computation_ 135.2-3 Elsevier, 2003, pp. 471–490
  * [17] Nikolai N Bogolyubov and Yurij A Mitropolskij  ‘‘Asymptotic Methods in the Theory of Non-linear Oscillations’’  Moscow: Nauka, 1963
  * [18] Y A Mitropolsky and O V Lykova  ‘‘Integral Manifolds in Nonlinear Mechanics (in Russian)’’  Moscow: Nauka, 1973
  * [19] Eugene I Butikov  ‘‘Analytical expressions for stability regions in the Ince–Strutt diagram of Mathieu equation’’  In _American Journal of Physics_ 86.4 American Association of Physics Teachers, 2018, pp. 257–267
  * [20] Sergey Vladimirovich Bolotin and Valery Vasil’evich Kozlov  ‘‘Calculus of variations in the large, existence of trajectories in a domain with boundary, and Whitney’s inverted pendulum problem’’  In _Izvestiya: Mathematics_ 79.5 IOP Publishing, 2015, pp. 894