11institutetext: Department of Computer Science, Johns Hopkins University, USA
11email<EMAIL_ADDRESS>11email<EMAIL_ADDRESS>
22institutetext: Institute of Natural Sciences, Shanghai Jiao Tong University,
China
33institutetext: School of Mathematical and Physical Sciences, Macquarie
University, Australia

# ProteinEngine: Empower LLM with Domain Knowledge for Protein Engineering

Yiqing Shen 11 Outongyi Lv 22 Houying Zhu 33 Yu Guang Wang(✉) 22

###### Abstract

Large language models (LLMs) have garnered considerable attention for their
proficiency in tackling intricate tasks, particularly leveraging their
capacities for zero-shot and in-context learning. However, their utility has
been predominantly restricted to general tasks due to an absence of domain-
specific knowledge. This constraint becomes particularly pertinent in the
realm of protein engineering, where specialized expertise is required for
tasks such as protein function prediction, protein evolution analysis, and
protein design, with a level of specialization that existing LLMs cannot
furnish. In response to this challenge, we introduce ProteinEngine, a human-
centered platform aimed at amplifying the capabilities of LLMs in protein
engineering by seamlessly integrating a comprehensive range of relevant tools,
packages, and software via API calls. Uniquely, ProteinEngine assigns three
distinct roles to LLMs, facilitating efficient task delegation, specialized
task resolution, and effective communication of results. This design fosters
high extensibility and promotes the smooth incorporation of new algorithms,
models, and features for future development. Extensive user studies, involving
participants from both the AI and protein engineering communities across
academia and industry, consistently validate the superiority of ProteinEngine
in augmenting the reliability and precision of deep learning in protein
engineering tasks. Consequently, our findings highlight the potential of
ProteinEngine to bride the disconnected tools for future research in the
protein engineering domain.

###### Keywords:

Deep Learning Large Language Model Protein Design AI for Protein Design.

##  1 Introduction

Large language models (LLMs) have achieved remarkable successes in solving
complex tasks, showcasing their zero-shot learning capabilities [18]. However,
the effectiveness of these models tends to plateau when faced with more
specialized tasks due to their inability to access domain-specific knowledge
or utilize specialized tools tailored for certain applications. This
limitation becomes glaringly apparent in the context of protein engineering
tasks. Although LLMs have been explored for specific tasks within the protein
engineering domain such as protein structure prediction [11], protein
evolution analysis, or de novo protein design [6, 13], their application often
demands significant alterations to the model architecture. These
modifications, coupled with the need for learning from scratch using domain-
specific datasets, present two-fold challenges. First, this approach
underutilizes the capabilities of well-trained foundation models, given their
robust pre-existing knowledge base. Second, the process of specialization
often leads to the loss of the model’s conversational abilities, a key feature
that makes LLMs versatile and user-friendly.

The introduction of in-context learning capabilities in LLMs [4] has ushered
potential solutions to the enduring challenge of domain knowledge scarcity
[4]. In this new paradigm, the LLM acts as a centralized, cognitive-like
system, which is capable of addressing domain-specific tasks by invoking
relevant Application Programming Interfaces (APIs) or systems to bride the
knowledge gap within the specific domain. However, despite these advancements,
most existing solutions tend to restrict the tools they incorporate to
commonly used APIs such as calculators, calendars, and web searches, or
relatively simple AI models like text-to-image generation models. In the
context of protein engineering, both the task formulation and the APIs
involved, as well as the AI models, manifest greater complexity. They are
characterized by a diverse modality of input and a larger set of arguments,
accentuating the need for more sophisticated and flexible systems.
Consequently, further research and development are required to fully leverage
the potential of LLMs in complex domains such as protein engineering. To
narrow the gap, we present ProteinEngine, a novel LLM system for protein
engineering. The major contributions are three-fold:

  1. (1)

We introduce ProteinEngine, a human-centered platform to augment the
capabilities of LLMs in tackling protein engineering tasks. This is achieved
by seamlessly integrating a comprehensive array of tools, packages, and
software relevant to protein engineering, all accessible through APIs.

  2. (2)

We propose a role-playing framework, comprising AI Project Manager (AI-PM), AI
Domain Expert (AI-DE), and AI Presenter (AI-Pr) modules, which facilitates
efficient task delegation, promotes interdisciplinary integration, ensures
dynamic adaptability, and enables effective communication of results
respectively. This design principle provides substantial flexibility, allowing
for easy integration and extensibility with emerging AI protein-design models.

  3. (3)

Through comprehensive user studies, we demonstrate the superior performance of
ProteinEngine in not only enhancing the usability of currently disconnected
protein engineering tools but also reducing the workload and learning
difficulty across users with different backgrounds.

Figure 1:  The overall framework of the proposed ProteinEngine, which
incorporates three distinct roles, each assigned to a separate LLM.

##  2 ProteinEngine: LLM for Protein Engineering

### Method Overview

Our proposed ProteinEngine is a human-centered system designed to augment the
capabilities of existing LLMs to address a broader spectrum of intricate
protein engineering tasks. Specifically, we assign three critical roles to the
LLMs: the AI Project Manager (AI-PM), the AI Domain Expert (AI-DE), and the AI
Presenter (AI-Pr). The overall pipeline illustrating these roles and their
interactions is depicted in Figure 1. The AI-PM functions as the primary
coordinator, interpreting user input expressed in natural language and
ensuring all necessary inputs, configurations, arguments, and conditions are
correctly provided. Subsequently, it breaks down a given complex task into
smaller, more manageable sub-tasks, delegating them to the appropriate AI-DEs.
Then, we employ multiple AI-DEs within the platform to address the wide
variety of challenges inherent in protein engineering and to facilitate future
expansion. Each AI-DE specializes in a particular domain or category of tasks,
ensuring a comprehensive coverage of the diverse aspects of protein
engineering. During the inference stage, the AI-PM selects and assigns a
subset of all AI-DEs to execute relevant APIs based on the nature of the task
at hand. Lastly, the AI-Pr is tasked with presenting the results, either
unimodal or multimodal, generated by the AI-DEs to the user.

### AI Project Manager

The LLM performing as the AI-PM acts as the primary interface, bridging the
gap between the user and the underlying protein engineering tools within the
ProteinEngine platform. In its core role, the AI-PM is tasked with
interpreting user input presented in natural language, discerning the context,
and identifying the necessary tasks to be performed. Beyond this, the AI-PM
ensures that all required inputs, arguments, configurations, and conditions
are correctly provided. To accurately parse and deconstruct the user’s query,
the AI-PM uses in-context learning, a more efficient alternative to the
computationally demanding process of LLM fine-tuning. By systematically
decomposing complex tasks into smaller, manageable sub-tasks, the AI-PM
ensures a thorough understanding of the user’s requirements. Once the sub-
tasks are defined, the AI-PM delegates them to the appropriate AI-DEs, taking
into account their respective areas of specialization.

### AI Domain Expert

The AI-DE in the ProteinEngine is specifically designed to manage a distinct
category of tasks pertaining to protein engineering. A team of multiple AI-DEs
is assembled, with each expert equipped with the necessary domain-specific
knowledge and tools to execute its designated tasks. To ensure the AI-DEs
perform efficiently and adaptively, we have implemented a novel self-feedback
communication loop mechanism between the AI-DE and the AI-PM. This autonomous
mechanism, which operates without the need for human intervention, enables AI-
DEs to progressively refine their understanding of new challenges that may
arise during the execution process, and to seek assistance from their fellow
AI-DEs, if required. As a result, AI-DEs can dynamically adjust and respond to
the evolving demands of the tasks, thereby maintaining a high degree of
accuracy and effectiveness.

### AI Presenter

AI-Pr aggregates and presents the results generated by the AI-DEs in a clear,
concise, and user-friendly manner, ensuring the user can easily interpret and
utilize the generated insights, fostering a deeper understanding. To
effectively communicate the results to the user, the AI-Pr is capable of
visualizing multimodal data, which includes, but is not limited to, images and
textual data. This presentation is tailored to cater to different user
preferences, and it enhances the comprehensibility of the data, enabling users
to quickly grasp the key insights and outcomes delivered by the AI-DEs.

Table 1: The involved AI models and APIs for protein engineering in the
proposed ProteinEngine.

API | Functionality | Description | Input | Output
---|---|---|---|---
AlphaFold 2 [9] | protein folding |  single-chain structure prediction with MSA | protein sequence |  atom-level 3D coordinates; residue-level pLDDT
AlphaFold-Multimer [5] | protein folding |  multi-chain structure prediction with MSA | protein sequence |  atom-level 3D coordinates; residue-level pLDDT
ESMFold [11] | protein folding |  MLM-based structure prediction without MSA | protein sequence |  atom-level 3D coordinates; residue-level pLDDT
MSA Transformer [14] | protein folding |  single-chain structure prediction with MSA | multiple sequence alignment |  atom-level 3D coordinates; in .pdb format
ESM-IF1 [11] | inverse folding |  single-site mutation Transformer-based | protein sequence | de novo protein sequence
LGN [19] | variant effect prediction |  deep mutation GNN based denoising | protein graph | de novo protein sequence
Equidock [7] | protein-target interaction | rigid-body docking |  two protein structures in .pdb format | binding affinity score
EquiBind [15] | protein-target interaction | rigid-body docking |  protein-ligand structure in .pdb format |  protein-ligand interaction sites binding affinity score
DiffDock [3] | protein-target interaction | rigid-body docking |  antibody-antigen structures in .pdb format | bound structure of complex
Diffab [12] | protein target interaction | antibody-antigen interaction |  antibody-antibody structures in .pdb format |  bound structure of complex binding affinity, epitope mapping
ProtENN [1] | sequence generation |  language-based model | protein sequence and structure | de novo protein sequence with function
Progen [13] | sequence generation |  language-based model | protein sequence | de novo protein sequence with function
Grade-IF [16] | sequence generation |  graph-based model | protein sequence | de novo protein sequence
GearNet [17] | property prediction |  latent representation of protein structure | protein graph | function or structure label
DeepSol [10] | property prediction |  solubility prediction | protein sequence | solubility
PyMOL | protein visualization |  visualize 3D conformation for a given protein molecule |  .pdb document | 3D visualization of the protein
VMD | protein visualization |  visualize 3D conformation for a given protein molecule |  .pdb document | 3D visualization of the protein
BioMedLM | biomedical domain Q&A |  trained on biomedical literature and clinical notes |  natural language on biomedicine or healthcare |  answer questions as a specialist in the field

### APIs for Protein Engineering

The APIs integrated within our system, along with their corresponding task
category taxonomy, are delineated in Table 1, which covers most of the task
scenarios in protein design. Each category is mapped to a specific AI-DE. To
provide a tangible understanding of the ProteinEngine in action, we illustrate
its use through distinct case examples in Figure 2, where different AI-DEs are
involved in each case. We illustrate typical examples of SOTA AI models in
protein design which have been used in our ProteinEngine platform.

Figure 2:  Three representative use case examples of the ProteinEngine in user
mode, where only the absence of mandatory parameters will be requested to the
user.

##  3 User Study

Figure 3:  The overall flowchart of the user study. This includes preparation
(participant recruitment and briefing), user participation operation
(familiarization with the technology and random assignment to conditions),
data collection (sequential tasks under different conditions with intermittent
feedback). Each stage of the process is color-coded for ease of understanding.

### Hypothesis Formulation and Testing

To evaluate the effectiveness of ProteinEngine, we conducted a user study
focused on gauging its proficiency as an intuitive, human-centered system for
protein engineering tasks. We employed hypothesis testing to quantitatively
compare the performance of ProteinEngine against a baseline condition,
focusing on task completion time, number of attempts, system usability, and
the perceived workload. In the baseline condition, participants employed
traditional tools and methods, independent of ProteinEngine, such as executing
Python scripts directly. Therefore, our null hypotheses were formulated as
follows:

  1. (H1)

ProteinEngine does not reduce the time required for successful identification
and execution of protein engineering models against the baseline.

  2. (H2)

ProteinEngine does not improve the accuracy in identifying and executing
models for protein engineering tasks against ProteinEngine does improve the
accuracy against the baseline.

  3. (H3)

ProteinEngine does not enhance the overall system usability for model
identification and execution within protein engineering tasks against the
baseline.

  4. (H4)

ProteinEngine does not decrease the workload required for the completion of
protein engineering model identification and execution against the baseline.

Based on the hypotheses, we used a single-sided t-test to assess statistical
significance.

### Dependent Variables

Statistical tests on these hypotheses involve collecting data on the dependent
variables from the user study, specifically the task completion time, number
of attempts, usability score, and workload index, which are defined as
follows.

  * •

Task Completion Time: This objective, continuous variable measures the total
time each participant takes to successfully complete a task under each
condition.

  * •

Number of Attempts: Another objective, continuous variable records the total
attempts a participant takes to successfully complete each task under each
condition. This variable is indicative of the accuracy of user actions.

  * •

Usability Score: This subjective, continuous variable is derived from the
System Usability Scale (SUS) questionnaire [2]. The score reflects the
perceived usability of the system.

  * •

Workload Index: This subjective, continuous variable, sourced from the NASA
Task Load Index (NASA TLX) questionnaire [8], assesses the perceived mental
workload across six dimensions: mental demand, physical demand, temporal
demand, effort, performance, and frustration level.

### Independent Variables

The primary independent variable is the Condition under which participants
perform the protein engineering tasks, either the baseline or the
ProteinEngine. Additionally, we consider potential confounding independent
variables that could influence our study outcomes, including:

  * •

Participant Background: These categorical variables encapsulate information
about each participant’s professional role and affiliations (academia or
industry). This information could offer insights into a user’s likely
background knowledge and potential biases or preferences when using the
interface.

  * •

Familiarity with Technology: These numerical variables represent the degree of
each participant’s familiarity with protein engineering tasks, Python
programming language, AI models, and the intersection of these areas. The
level of familiarity could potentially influence the ease with which
participants adapt to the ProteinEngine, and thus might impact the
measurements of variables like task completion time, number of attempts, and
perceived usability and workload.

### User Study Design

The user study workflow, shown in Figure 3, consisted of several steps. First,
participants received an introductory tutorial that provided information about
the study’s background, motivation, and procedures. Next, participants
completed a preliminary questionnaire that assessed their background knowledge
and familiarity with AI, protein engineering, and their interdisciplinary
overlap. The study followed a between-subjects design, comparing participants
exposed to two conditions: a baseline condition and the ProteinEngine
condition. To mitigate learning effects, the order in which participants
encountered these conditions was randomized. To control for potential effects
stemming from participants’ background and familiarity with technology, all
participants completed the same set of activities under both conditions
(baseline and ProteinEngine). This approach provided paired data for analysis.
Participants were assigned a series of six distinct protein engineering tasks
under each experimental condition. These tasks included protein folding,
inverse protein folding, and protein mutation prediction. During the task
completion process, we carefully recorded the total time taken for each task
and the number of attempts required for successful execution. After completing
the tasks in either the baseline or ProteinEngine condition, participants were
asked to fill out a questionnaire. This questionnaire aimed to assess their
subjective impressions of the system’s usability and their perceived workload
during task completion. To ensure a valid comparison of user experiences,
identical questionnaires were administered. Data was collected using Google
Sheets, with all questions being mandatory to prevent missing data. Incomplete
data from participants who withdrew or failed to complete tasks were excluded
from the analysis.

### Participants Recruitment

We strategically planned the recruitment of volunteer participants to
encompass a wide range of potential users. This design aimed to test the
versatility and broad applicability of our proposed method across various user
groups. Our participant cohort consisted of volunteers from both the AI and
biological communities, spanning both academic and industrial fields. This
diverse group, including students, AI researchers, lab technicians, and
biologists, allows for a comprehensive evaluation of ProteinEngine’s
functionality across multiple user profiles.

### Implementations

For the LLM in ProteinEngine, we chose the gpt-3.5-turbo. As the most advanced
model in the GPT-3.5 series, this version provides robust capabilities and
superior performance suitable for our application.

### Collected Data

Figure 4:  Boxplots for the four variables.

The boxplots illustrating the distribution of our four key variables, namely
Task Completion Time, Number of Attempts, Usability Score, and Workload Index,
can be found in Figure 4. Each box plot provides a visual summary of the
minimum, first quartile (Q1), median (Q2), third quartile (Q3), and maximum
values for these variables. The box represents the interquartile range (IQR)
from Q1 to Q3, the line inside the box denotes the median, and the whiskers
extend to show the range of the data within 1.5 times the IQR. Observations
beyond this range are considered outliers and are represented as individual
points. On average it take 36.40 minutes to complete the user study for each
participant.

### Hypothesis Testing

We applied the hypothesis testing to ascertain whether the performance
differences observed between the two conditions i.e., baseline and
ProteinEngine were statistically significant. The differences under
consideration, represented as $d_{i}$, were computed by subtracting the
ProteinEngine measurements from the baseline measurements. Under the null
hypothesis, where both platforms have an equivalent effect, these differences
should follow a distribution centered around zero, i.e., $\mu_{d}=0$. Our
final dataset for hypothesis testing comprised $n=26$ samples. We formulated
our null and alternative hypotheses as follows:
$H_{0}:\mu_{d}=0\quad\texttt{against}\quad H_{1}:\mu_{d}>0$. This holds for
testing hypotheses H1, H2, and H4. For testing H3, the alternative hypothesis
is $\mu_{d}<0$. Let $\overline{d}$ and $s_{d}$ represent the sample mean and
sample standard deviation of the observed differences, respectively. Given
these parameters, the sampling distribution for the test statistic follows a
$t$ distribution with degrees of freedom $n-1$. This means that, under null
hypothesis $H_{0}$, $\tau=\frac{\overline{d}}{s_{d}/\sqrt{n}}\sim t_{n-1}$.
Table 2 compiles our hypothesis testing results. Of note, no significant
difference was found for the workload measure between the two conditions.

Table 2: Null hypothesis testing results. Hypothesis | Observed test statistics | P-value | Reject null hypothesis
---|---|---|---
H1 | $7.7012$ | $2.335\times 10^{-8}$ | Yes
H2 | $3.1944$ | $1.884\times 10^{-3}$ | Yes
H3 | $-2.4162$ | $1.166\times 10^{-2}$ | Yes
H4 | $0.74029$ | $2.330\times 10^{-1}$ | No

### Results

Our study encompassed a total of 45 participants. For hypothesis testing, we
employed the paired two-sample t-test (effectively a one-sample, one-sided
t-test on the difference) at a 5% significance level for the four variables
under two conditions, baseline and ProteinEngine. The results allowed us to
reject three of the null hypotheses, thereby highlighting the superior
performance of ProteinEngine in facilitating protein engineering tasks.

##  4 Conclusion

In this work, we presented ProteinEngine, a groundbreaking platform that
amplifies the capabilities of LLMs in the realm of protein engineering. This
platform’s human-centered design greatly eases the learning curve
traditionally associated with specialized tools, thereby making protein
engineering tasks more accessible. By integrating advanced LLMs with domain-
specific expertise, ProteinEngine marks a significant leap forward in the
application of AI to protein engineering, showing great potential to
accelerate scientific discoveries and spur innovation. As we continue to
develop and refine ProteinEngine, it is critical to emphasize the importance
of responsible use and rigorous validation. Therefore, the development of
comprehensive ethical guidelines and robust validation protocols is a key
direction for future work.

## References

  * [1] Maxwell L Bileschi et al.  Using deep learning to annotate the protein universe.  Nature Biotechnology, 40(6):932–937, 2022.
  * [2] John Brooke et al.  Sus-a quick and dirty usability scale.  Usability evaluation in industry, 189(194):4–7, 1996.
  * [3] Gabriele Corso et al.  DiffDock: diffusion steps, twists, and turns for molecular docking.  In The Eleventh International Conference on Learning Representations, 2023.
  * [4] Qingxiu Dong et al.  A survey for in-context learning.  arXiv preprint arXiv:2301.00234, 2022.
  * [5] Richard Evans et al.  Protein complex prediction with AlphaFold-Multimer.  bioRxiv, pages 2021–10, 2022.
  * [6] Noelia Ferruz et al.  Protgpt2 is a deep unsupervised language model for protein design.  Nature communications, 13(1):4348, 2022.
  * [7] Octavian-Eugen Ganea et al.  Independent SE(3)-equivariant models for end-to-end rigid protein docking.  In International Conference on Learning Representations, 2021.
  * [8] Sandra G Hart et al.  Development of nasa-tlx (task load index): Results of empirical and theoretical research.  In Advances in psychology, volume 52, pages 139–183. Elsevier, 1988.
  * [9] John Jumper et al.  Highly accurate protein structure prediction with alphafold.  Nature, 596(7873):583–589, 2021.
  * [10] Sameer Khurana et al.  Deepsol: a deep learning framework for sequence-based protein solubility prediction.  Bioinformatics, 34(15):2605–2613, 2018.
  * [11] Zeming Lin et al.  Evolutionary-scale prediction of atomic-level protein structure with a language model.  Science, 379(6637):1123–1130, 2023.
  * [12] Shitong Luo et al.  Antigen-specific antibody design and optimization with diffusion-based generative models for protein structures.  In Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho, editors, Advances in Neural Information Processing Systems, 2022.
  * [13] Ali Madani et al.  Large language models generate functional protein sequences across diverse families.  Nature Biotechnology, pages 1–8, 2023.
  * [14] Roshan M Rao et al.  Msa transformer.  In Marina Meila and Tong Zhang, editors, Proceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pages 8844–8856. PMLR, 18–24 Jul 2021.
  * [15] Hannes Stärk et al.  Equibind: geometric deep learning for drug binding structure prediction.  In International Conference on Machine Learning, pages 20503–20521. PMLR, 2022.
  * [16] Kai Yi et al.  Graph denoising diffusion for inverse protein folding.  Advances in Neural Information Processing Systems, 36, 2024.
  * [17] Zuobai Zhang et al.  Protein representation learning by geometric structure pretraining.  In The Eleventh International Conference on Learning Representations, 2023.
  * [18] Wayne Xin Zhao et al.  A survey of large language models.  arXiv preprint arXiv:2303.18223, 2023.
  * [19] Bingxin Zhou et al.  Accurate and definite mutational effect prediction with lightweight equivariant graph neural networks.  arXiv preprint arXiv:2304.08299, 2023.


# Identification of Secondary Resonances of Nonlinear Systems using Phase-
Locked Loop Testing

Tong Zhou<EMAIL_ADDRESS>Gaëtan Kerschen  Space Structures and Systems
Laboratory, Aerospace and Mechanical Engineering Department, University of
Liege, Belgium

###### Abstract

One unique feature of nonlinear dynamical systems is the existence of
superharmonic and subharmonic resonances in addition to primary resonances. In
this study, an effective vibration testing methodology is introduced for the
experimental identification of these secondary resonances. The proposed method
relies on phase-locked loop control combined with adaptive filters for online
Fourier decomposition. To this end, the concept of a resonant phase lag is
exploited to define the target phase lag to be followed during the
experimental continuation process. The method is demonstrated using two
systems featuring cubic nonlinearities, namely a numerical Duffing oscillator
and a physical experiment comprising a clamped-clamped thin beam. The obtained
results highlight that the control scheme can accurately characterize
secondary resonances as well as track their backbone curves. A particularly
salient feature of the developed algorithm is that, starting from the rest
position, it facilitates an automatic and smooth dynamic state transfer toward
one point of a subharmonic isolated branch, hence, inducing branch switching.

###### keywords:

Phase-locked loop testing , Superharmonics , Subharmonics , Mode
identification

##  1 Introduction

Experimental modal analysis is the most common approach for identifying the
dynamical characteristics of a mechanical system. Specifically, the objective
is to determine the modal parameters, namely the vibration modes, resonance
frequencies and damping ratios, from the analysis of free or forced
vibrations. At the root of experimental modal analysis is the assumption of
linear dynamical behavior. However, nonlinear phenomena are more and more
frequently encountered during vibration testing of modern engineering
structures [1]. Typical examples include commercial aircraft [2], satellites
[3], turbomachinery [4] and microelectromechanical systems (MEMS) [5]. This is
why there is currently a burgeoning interest in developing methodologies that
can tackle nonlinear behaviors.

Nonlinear systems are known for exhibiting rich and sophisticated dynamical
phenomena which include the dependence of resonance frequency on motion
amplitude and the resulting hardening and softening phenomena, co-existence of
multiple stable vibration states and jumps between different states due to
external perturbations, existence of unstable branches in the nonlinear
frequency response curves (NFRCs) that are interconnected with stable portions
at bifurcation points, activation of multiharmonic responses by a mono-
harmonic excitation, appearance of secondary resonances, namely superharmonic
and subharmonic resonances, in addition to the primary resonance, emergence of
isolated resonance curves detached from the main branch, occurrence of
internal resonances and exchange of energy between different modes, etc [6].
These dynamical phenomena pose great challenges for the development of robust
vibration testing schemes.

In this context, nonlinear normal modes (NNMs), defined either as periodic
motions of an undamped and unforced mechanical system [7, 8, 9], invariant
manifolds in phase space [10, 11] or spectral submanifolds [12], have become a
central concept in the structural dynamics community. NNMs offer valuable
insights into the dynamics of a nonlinear system near resonance [13]. For
instance, the NNM theory proved useful for the analysis of slender structures
featuring geometrical nonlinearity as, e.g., in [14]. In analogy to the
frequency response function of a linear system, the NFRCs are obtained by
computing branches of periodic responses of a damped, harmonically-forced
system. The backbone curves trace the locus of the resonance peaks,
establishing a direct link between NNMs and NFRCs.

One of the first attempts to identify NNMs experimentally exploited a
nonlinear extension of the phase lag quadrature criterion [15, 16] which, for
a linear system, states that a normal mode is excited when the displacement
and harmonic forcing are in quadrature. This is at the root of nonlinear phase
resonance testing. Following NNM appropriation, their amplitude-dependent
properties can, in turn, be extracted during free decay oscillation according
to the invariance property of NNMs. More recently, phase-locked loop testing
(PLLT) was introduced to robustify NNM identification [17, 18, 19, 20, 21].
PLLT makes use of feedback-controlled excitation to impose a user-specified
phase lag between the applied excitation and the measured response. To this
end, the excitation frequency is automatically adjusted until the control
objective is met. An NFRC can be obtained by sequentially increasing or
decreasing the phase lag at constant forcing amplitude to sweep the resonance
peak whereas a backbone curve can be measured by sequentially increasing the
forcing amplitude while ensuring phase quadrature conditions.

The vast majority of the studies in the technical literature dealing with
experimental modal analysis of nonlinear structures focuses on the
characterization of primary resonances. However, a fundamental harmonic
excitation can activate fractional and higher-order harmonic responses, which,
in turn, can trigger the excitation of secondary resonances. Even if the
superharmonic and subharmonic resonances of nonlinear systems have been widely
studied analytically and numerically, see, e.g., [6, 22, 23], their
experimental identification remains a field to be explored. In this context,
we note the pioneering works by Chauhan et al. [24] and Yamaki et al. [25,
26]. Experimental studies on other physical systems also reported the
existence of subharmonic resonance phenomena, such as acoustic emission from
gas bubbles [27], nonlinear electrical LC circuits [28], a moored ocean
structural system [29], shear layer flow under excitation [30]. Reference [31]
considers the monitoring of subharmonic response components for detecting
bolted joint loosening.

The objective of this study is to develop a vibration testing scheme based on
PLL feedback control and the concept of a resonant phase lag [32] for the
identification of the secondary resonances of a nonlinear mechanical system.
To this end, PLL is coupled to an adaptive digital filter, a powerful building
block which was first incorporated in the control loop in [33]. The adaptive
filter allows the analyst to perform online Fourier decomposition and, hence,
estimate the phase lag of a target harmonic of the response with respect to
the harmonic forcing. The resonant phase lag of the considered secondary
resonance is then exploited for designing the appropriate control strategy to
identify its NFRC and backbone curve.

The paper is organized as follows. The basic principle of the proposed method
is introduced in Section 2. In Section 3, the control scheme is tested on a
Duffing oscillator in a virtual experiment. The performance of PLLT for the
characterization of folded superharmonic responses and subharmonic isolas is
investigated. In Section 4, a slender beam exhibiting geometrical nonlinearity
is employed as the test system in a physical experiment. Conclusions are drawn
in the final section.

##  2 Phase-locked loop testing with adaptive filters

The basic principle of PLLT is to make use of the unique phase variation
property of a nonlinear resonance to unfold the NFRCs. This is exemplified for
the nonlinear primary resonance of a hardening Duffing oscillator excited by a
sine forcing function in Figure 1. It is seen that the phase lag between the
fundamental harmonic of the displacement and the forcing reduces monotonically
from 0 to $-\pi$ when increasing the forcing frequency. Thus, there is a
unique forcing frequency corresponding to a specific value of the phase lag.
In addition, the backbone curve describing the locus of the resonance peaks
can be tracked by locking the phase angle at a value of $-\pi/2$ while varying
the force amplitude.

The PLL scheme utilized in this work requires no prior knowledge of the tested
system and can be easily implemented in a physical experiment. Another benefit
is that it can stabilize unstable branches through the action of the
controller [19, 34]. As schematized in Figure 2, it comprises there important
building blocks, namely a phase detector, a proportional-integral (PI)
controller and a voltage-controlled oscillator (VCO). A closed-loop system is
built after the inclusion of the mechanical system to be tested. The system
state is prescribed through the imposition of a desired phase lag between the
harmonic excitation and a specific harmonic of the response. To this end, the
forcing frequency is automatically adjusted by the controller. The harmonic
excitation signal is generated by the VCO module based on the fact that the
circular frequency is a phase change per unit time. The external forcing is
thus expressed as

$\displaystyle
f(t)=F\sin{\left[\int_{0}^{t}\omega(\tau)\mathrm{d}\tau\right]}.$ (1)

The instantaneous frequency is determined by the difference between the
assigned reference phase $\Phi_{\mathrm{ref}}$ and the phase lag of the
$\kappa$-th harmonic $\Phi_{\kappa,\upsilon}$ obtained by the phase detector
and the settled PI controller

$\displaystyle\omega(t)=\omega_{\mathrm{o}}+K_{P}\left[\Phi_{\kappa,\upsilon}(t)-\Phi_{\mathrm{ref}}\right]+K_{I}\int_{0}^{t}\left[\Phi_{\kappa,\upsilon}(\tau)-\Phi_{\mathrm{ref}}\right]\mathrm{d}\tau.$
(2)

Here, $\omega_{\mathrm{o}}$ is the open-loop forcing frequency, often set to
be close to the target resonance. $K_{P}$ and $K_{I}$ are the proportional and
integral gains of the controller, respectively. They are employed to achieve
fast convergence to the target phase angle and eliminate the steady-state
errors in the feedback loop. Alternatively, Eq. (2) can be rewritten into a
set of differential equations [19]. A differential gain can be introduced to
further accelerate the convergence, but it is not considered in this work for
simplicity.

Figure 1: NFRCs of a hardening Duffing oscillator computed at different
forcing amplitudes by the harmonic balance method. The unstable responses are
indicated by dashed lines. (a) Amplitude and (b) phase lag. Figure 2: PLL
scheme with adaptive filters.

The measured periodic response is written in terms of the truncated Fourier
expansion

$\displaystyle
x(t)=\dfrac{c_{0,\upsilon}}{\sqrt{2}}+\sum_{\kappa=1}^{N}\left[s_{\kappa,\upsilon}\sin{\left(\dfrac{\kappa\omega}{\upsilon}t\right)}+c_{\kappa,\upsilon}\cos{\left(\dfrac{\kappa\omega}{\upsilon}t\right)}\right].$
(3)

where $\kappa$ is the number of harmonics and $\upsilon$ is a positive integer
introduced to account for subharmonics. The fundamental harmonic excitation
can activate multiharmonic responses which can consequently trigger the
excitation of secondary resonances, namely $\kappa:1$ superharmonic
resonances, $1:\upsilon$ subharmonic resonances or $k:\upsilon$
ultrasubharmonic resonances [32], when the frequency of the $\kappa$-th
harmonic $\omega_{\kappa,\upsilon}=\kappa\omega/\upsilon$ corresponds to that
of 1:1 primary resonance. The expansion (3) can be recast into a more compact
form

$\displaystyle x(t)=\textbf{Q}(t)\textbf{z},$ (4)

where $\textbf{Q}(t)$ is a row vector collecting the time-dependent
trigonometric terms

$\displaystyle\textbf{Q}(t)=\left[\dfrac{1}{\sqrt{2}}\;\sin{\left(\dfrac{\omega}{\upsilon}t\right)}\;\cos{\left(\dfrac{\omega}{\upsilon}t\right)}\;...\;\sin{\left(\dfrac{N\omega}{\upsilon}t\right)}\;\cos{\left(\dfrac{N\omega}{\upsilon}t\right)}\right],$
(5)

and z is a vector containing the associated Fourier coefficients

$\displaystyle\textbf{z}=\left[c_{0,\upsilon}\;s_{1,\upsilon}\;c_{1,\upsilon}\;...\;s_{N,\upsilon}\;c_{N,\upsilon}\right]^{\mathrm{T}},$
(6)

In this work, the amplitude and phase lag of the $\kappa$-th harmonic
component are determined as

$\displaystyle
x_{\kappa,\upsilon}(t)=A_{\kappa,\upsilon}^{x}\sin{\left(\dfrac{\kappa\omega}{\upsilon}t+\Phi_{\kappa,\upsilon}^{x}\right)},$
(7)

where

$\displaystyle
A_{\kappa,\upsilon}^{x}=\sqrt{c_{\kappa,\upsilon}^{2}+s_{\kappa,\upsilon}^{2}},$
(8)

and

$\displaystyle\Phi_{\kappa,\upsilon}^{x}=\mathrm{arctan2}(c_{\kappa,\upsilon},s_{\kappa,\upsilon}).$
(9)

To develop a control-based scheme that can effectively characterize secondary
resonances, the phase of the non-fundamental harmonic of interest, i.e.,
$\kappa\neq\upsilon$, must be accurately estimated. This is why we resort to
adaptive digital filters for performing online Fourier decomposition [35].
Indeed, reference [34] demonstrated that adaptive filtering presents distinct
advantages over the conventional synchronous demodulation techniques based on
low-pass filters. Adaptive filtering synthesizes a signal $\hat{x}(t)$ by
considering a time-dependent basis weighted by a set of coefficients, such
that it approximates the measured response $x(t)$

$\displaystyle\hat{x}(t)=\textbf{Q}(t)\hat{\textbf{z}}\approx x(t),$ (10)

The row vector $\textbf{Q}(t)$ containing the sine and cosine functions now
becomes the basis in adaptive filtering, and the weight coefficients
$\hat{\textbf{z}}$ are combined in the vector

$\displaystyle\hat{\textbf{z}}=\left[\hat{z}_{0}\;\hat{z}_{1s}\;\hat{z}_{1c}\;...\;\hat{z}_{Ns}\;\hat{z}_{Nc}\right]^{\mathrm{T}}.$
(11)

These weights can be interpreted as an approximation of the Fourier
coefficients with $\hat{\textbf{z}}\approx\textbf{z}$. They are updated online
at discrete time steps using the least mean squares (LMS) algorithm, which is
more cost-effective as compared to other alternative algorithms [33]. The
weights at the instant $t_{i+1}$ are updated according to the gradient of the
mean squared error at $t_{i}$

$\displaystyle\hat{\textbf{z}}(t_{i+1})=\hat{\textbf{z}}(t_{i})+\mu\epsilon(t_{i})\textbf{Q}^{\mathrm{T}}(t_{i}),$
(12)

where $\mu$ is the step size factor specifying the convergence and $\epsilon$
is the synthesis error estimated using

$\displaystyle\epsilon(t_{i})=x(t_{i})-\hat{x}(t_{i})=x(t_{i})-\textbf{Q}(t_{i})\hat{\textbf{z}}(t_{i}),$
(13)

The phase of the $\kappa$-th harmonic of the displacement can then be
extracted according to

$\displaystyle\hat{\Phi}_{\kappa,\upsilon}^{x}=\mathrm{arctan2}(\hat{z}_{\kappa
c},\hat{z}_{\kappa s}).$ (14)

The phase of the fundamental forcing component at $\kappa=\upsilon$ denoted
$\hat{\Phi}_{\upsilon,\upsilon}^{f}$ can be decomposed in a similar manner, so
that the phase lag of the $\kappa$-th harmonic of the displacement with
respect to the forcing can be estimated using the relation

$\displaystyle\hat{\Phi}_{\kappa,\upsilon}=\hat{\Phi}_{\kappa,\upsilon}^{x}-\hat{\Phi}_{\upsilon,\upsilon}^{f}\approx\Phi_{\kappa,\upsilon}.$
(15)

When it comes to practical implementation, the phase angles of the
trigonometric terms of the basis in (10) are estimated through the integration
of the instantaneous frequency with respect to time, to ensure a smooth phase
variation.

Eventually, the concept of a resonant phase lag introduced in [32] can be
exploited for the identification of the NFRC and backbone of the resonance of
interest. For instance, for a Duffing oscillator subjected to a sine forcing
function, it is shown that the points where the phase lag is $-\pi/2$
($\kappa$ and $\upsilon$ odd) or $-3\pi/4\upsilon$ ($\kappa$ or $\upsilon$
even) define the locus of phase resonance of the $\kappa:\upsilon$ resonance.
An important remark is that the phase lags used in this study are different
from those in [32] owing to difference in the sign conventions used to define
phase lag. The resonant phase lags of other oscillators are given in [36].

##  3 Numerical demonstration using a Duffing oscillator

A Duffing oscillator under sine forcing governed by the equation of motion

$\displaystyle m\ddot{x}(t)+c\dot{x}(t)+kx(t)+k_{nl}x^{3}(t)=F\sin{\omega t},$
(16)

is considered in this section. $m$, $c$, $k$ and $k_{nl}$ are the mass,
damping, linear and nonlinear stiffness coefficients, respectively. They are
chosen to be equal to $m=1$ kg, $c=0.001$ kg/s, $k=1$ N/m and $k_{nl}=1$ N/m3.
The natural frequency of the underlying linear system is
$\omega_{l}=\sqrt{k/m}$.

The dynamical model integrating PLL control and the Duffing oscillator was
implemented in the Matlab/Simulink environment. The fixed-step ordinary
differential equation solver was employed for carrying out the time-domain
simulations. The parameters of the PLL controller and of the adaptive filter
are tabulated in Table 1 where $f_{s}$ is the sampling frequency.

Target resonance |  $K_{P}$ [1/s] |  $K_{I}$ [1/s2] |  $f_{s}$ [1/s] |  $\mu$ [-]
---|---|---|---|---
1:1 | $1.0$ | $5\times 10^{-3}$ | $2\times 10^{2}$ | $10^{-4}$
3:1 and 2:1 | $0.1$ | $10^{-4}$ | $2\times 10^{2}$ | $10^{-4}$
1:3 and 1:2 | $2.0$ | $10^{-2}$ | $10^{3}$ | $10^{-4}$
Table 1: Parameters used in the virtual experiment.

###  3.1 Primary resonance

The NFRC and backbone of the primary resonance were identified using PLLT at
three different forcing amplitudes $F$, namely $0.00005$, $0.0001$ and
$0.00015$ N. A reference solution was also computed using the harmonic balance
method (HBM) [37] with 14 harmonics. The results are plotted in Figure 3.
Clearly, PLLT can trace out very accurately both the stable and unstable
branches of the NFRC. In addition, the backbone curve corresponding to
nonlinear phase resonance conditions can also be tracked accurately by locking
the phase lag at $-\pi/2$ and varying the forcing amplitude.

Figure 3: NFRCs ( ) and the backbone curve ( ) of the primary resonance
identified by PLLT. The reference solution is provided by HBM in gray; light
gray corresponds to unstable responses. (a) Amplitude and (b) phase lag.

Figure 4: NFRCs ( ) at $F=0.3$N and backbone curves ( ) of the 3:1 and 2:1
resonances identified by PLLT. The reference solution is provided by HBM in
gray; light gray corresponds to unstable responses. The fold and branch-point
bifurcations are marked with circles and crosses, respectively. (a) Amplitude
and (b) phase lag.

Figure 5: NFRCs ( ) at $F=0.5$N and backbone curves ( ) of the 3:1 and 2:1
resonances identified by PLLT. The reference solution is provided by HBM in
gray; light gray corresponds to unstable responses. The fold and branch-point
bifurcations are marked with circles and crosses, respectively. (a) Amplitude
and (b) phase lag.

###  3.2 Superharmonic resonances

As shown in [32] for a hardening Duffing oscillator and in Figures 4 and 5 of
this paper, the phase lag $\Phi_{\kappa,1}$ of the $\kappa$:1 superharmonic
resonances evolves in a similar manner compared to the phase lag of the
primary resonance. Specifically, the phase lag decreases monotonically when
increasing the forcing frequency, a property which is not fulfilled by the
fundamental harmonic $\Phi_{1,1}$ of the superharmonic resonance. Again, this
paves the way for a simple continuation scheme during which the phase lag is
sequentially decreased. Furthermore, around the superharmonic resonance, the
superharmonic component $x_{\kappa}$ makes a significant contribution to the
total response $x$, so that $\Phi_{\kappa,1}$ can be accurately extracted.

The identification of the 3:1 and 2:1 superharmonic resonances was achieved
using PLLT by letting the phase lags $\Phi_{3,1}$ and $\Phi_{2,1}$ vary
between 0 and $-\pi$ or between $-\pi/2$ and $-\pi$, respectively. The
corresponding phase resonances occur at $-\pi/2$ and $-3\pi/4$, respectively.
Figures 4 and 5 present the resulting NFRCs and backbone curves at $F=0.3$N
and $F=0.5$N, respectively. An excellent identification of the two
superharmonic resonances was obtained using PLLT. For the 2:1 resonance, a
second solution was also found to exist between 0 and $\pi/2$, as discussed in
[32]. Other odd and even superharmonic resonances can be characterized in a
similar manner.

Figure 6: Transition from the rest position to the 2:1 superharmonic
resonance: $F=0.5$ N, $\omega_{\mathrm{o}}=0.5$ and
$\Phi_{\mathrm{ref}}=-3\pi/4$. (a) Time evolution of the phase lag,
instantaneous frequency and displacement and (b) trajectory in the
displacement-frequency plane.

The calculation of the 3:1 resonance is facilitated by the fact that it
appears in the direct continuation of the main branch of the NFRC. This is not
the case for the 2:1 resonance which bifurcates out of the main branch. As a
result, HBM necessitates either a branch switching strategy or an initial
guess to calculate the 2:1 branch. Interestingly, this is not necessary for
the proposed PLLT scheme, because the 2:1 branch can be reached from the rest
position, i.e., starting with trivial initial conditions
$x_{0}=\dot{x}_{0}=0$. As illustrated in Figure 6, after a transient period
during which the phase lag converges toward $-3\pi/4$, the system’s motion
reaches the 2:1 branch. We note that $\omega_{0}$ shopuld be chosen close
enough from the frequency range in which the resonance exists.

Figure 7: Transition from the rest position to the upper resonance point of
the 1:3 subharmonic resonance: $F=0.5$ N, $\omega_{\mathrm{o}}=3$ and
$\Phi_{\mathrm{ref}}=-\pi/2$. (a-b) Time evolution of the phase lag,
instantaneous frequency and displacement: (a) complete time interval and (b)
close-up on the initial times; (c-d) Trajectory in the displacement-frequency
plane: (c) $K_{I}=2$ and $K_{P}=0.01$ and (d) $K_{I}=1$ $K_{P}=0.01$. The
yellow branch corresponds to the main branch of the NFRC. Figure 8:
Unsuccessful transition from the rest position to the 1:3 subharmonic
resonance: $F=0.1$N, $\omega_{\mathrm{o}}=3$ and $\Phi_{\mathrm{ref}}=-\pi/2$.
Figure 9: Normalized frequency responses $X/F$ of 1:3 resonances computed by
HBM.  Figure 10: Basins of attraction near the red point in Figure 9 for
$F=0.5$. Blue and yellow colors correspond to the periodic responses on the
main branch of the NFRC and on the 1:3 isola, respectively.

###  3.3 Subharmonic resonances

Unlike superharmonic resonances, subharmonic resonances take the form of
isolated branches of periodic solutions detached from the main branch [6].
This is why they are more difficult to calculate during numerical continuation
and can also be missed when considering classical sine sweep testing. These
resonances only exist beyond a specific threshold in forcing amplitude and
appear near multiples of the primary resonance frequency.

As discussed in [32] and shown in Figure 11(b), the phase lag
$\Phi_{1,\upsilon}$ of the 1:$\upsilon$ subharmonic resonances presents a more
complex topology owing to the closed nature of an isola. However, when either
the lower or the upper branch of the isola is considered, there is still a 1:1
relationship between the phase lag and the forcing frequency. For the 1:3 and
1:2 subharmonic resonances, the resonant phase lag is
$\Phi_{\mathrm{ref}}=-\pi/2$ and $\Phi_{\mathrm{ref}}=-3\pi/8$, respectively.
Unlike superharmonic resonances, each subharmonic branch possesses two points
featuring the resonant phase lag, i.e., one at the lowest extremity of the
isola and the other one at the highest extremity [32].

As for the 2:1 resonance, a very interesting feature of PLLT is that
subharmonic resonances can be reached without any guess or user interaction.
This is illustrated for the 1:3 subharmonic resonance in Figure 7 for
$F=0.5$N. PLLT starts from the rest position $x_{0}=\dot{x}_{0}=0$ and imposes
the resonant phase lag $\Phi_{\mathrm{ref}}=-\pi/2$. A so-called state
transfer process to the upper resonance point of the 1:3 isola then occurs.
There are mainly three stages in this transient process:

  * i

An initial overshooting, a typical phenomenon when using a PID controller,
results in fast oscillations of the instantaneous frequency around $\omega=3$.
This overshooting depends on the controller parameters $K_{P}$ and $K_{I}$.
For instance, Figure 7(d) shows that a smaller $K_{P}$ results in less severe
oscillations. These oscillations provide the trigger necessary to the jump
onto the isola. Once the oscillations are stabilized, the system state is
attracted in the vicinity of the low-amplitude resonance point of the isola.

  * ii

A smooth transfer along the subharmonic branch toward the high-amplitude
resonant point takes place. The fact that the PLL controller has a general
tendency to reach the upper resonance point can be explained by the fact that
there is an extremely rapid change in the phase lag around the lower amplitude
resonance point (see around $\omega\approx 3$ in Figure 11), which complicates
convergence toward this point.

  * iii

When the measured phase lag converges to the reference phase lag, the upper
resonance point of the subharmonic resonance is reached.

Such an automatic branch switching process could be achieved at other force
amplitudes or with reference phase lags different from the resonant phase lag.
However, it also happened that the transition to the isola could not be
achieved. This is the case for $F=0.1$N in Figure 8 where the forcing
frequency of the closed-loop system and the phase lag undergo oscillations
that cannot be stabilized. Figure 9 which displays the response amplitude
normalized by the force amplitude $X/F$ around $\omega=3$ reveals that the
distance between the main branch of the NFRC and the 1:3 subharmonic isola
increases when the forcing level decreases, hence complicating the transfer
process. Further confirmation is given by comparing the basins of attraction
of the periodic solution on the main branch and of the periodic solution on
the isola for $F=0.5$N in Figure 10. Clearly, the subharmonic response
possesses a large basin of attraction; its area represents $64.9\%$ of the
total area. This ratio can increase up to $84.4\%$ when $F=0.7$N whereas it is
close to $0\%$ when $F=0.1$N.

Figure 11: Backbone curve ( ) of the 1:3 subharmonic resonance identified by
PLLT. The starting point is the upper resonance point at $F=0.8$N indicated by
a star marker. The subharmonic branches (gray) are computed by HBM. The main
branch of the NFRC at $F=0.5$ N is represented in yellow. (a) Amplitude and
(b) phase lag.

Once the transfer process to the upper point of the 1:3 subharmonic isolas has
been achieved, the backbone can be identified using PLLT with
$\Phi_{\mathrm{ref}}=-\pi/2$. Figure 11 shows the outcome of PLLT when the
identification is initiated for $F=0.8$N and the forcing level is subsequently
decreased.

The identification of the 1:2 subharmonic resonance in Figure 13 is achieved
in a similar manner. The transfer process from the rest position to the 1:2
isola is represented for $F=0.5$N in Figure 12. In the present case, the
system is attracted to a point relatively far away from the lower resonance
point and then travels along the isola up to the upper resonance point. The
backbone curve can be subsequently identified by setting
$\Phi_{\mathrm{ref}}=-3\pi/8$ and decreasing the forcing amplitude $F$, as
displayed for $F=0.8$N in Figure 13.

Figure 12: Transition from the rest position to the upper resonance point of
the 1:2 subharmonic resonance: $F=0.5$ N, $\omega_{\mathrm{o}}=2$ and
$\Phi_{\mathrm{ref}}=-3\pi/8$. (a) Time evolution of the phase lag,
instantaneous frequency and displacement and (b) trajectory in the
displacement-frequency plane. The yellow branch corresponds to the main branch
of the NFRC.

Figure 13: Backbone curve ( ) of the 1:2 subharmonic resonance identified by
PLLT. The starting point is the upper resonance point at $F=0.8$N indicated by
a star marker. The subharmonic branches (gray) are computed by HBM. The yellow
branch corresponds to the main branch of the NFRC at $F=0.5$N. (a) Amplitude
and (b) phase lag.

Large portions of the NFRC of the 1:3 and 1:2 subharmonic resonances can also
be identified by modifying $\Phi_{\mathrm{ref}}$ from the upper resonance
point while keeping $F$ constant. The upper and lower branches are obtained by
either increasing or decreasing the reference phase lag, respectively. The
results for $F=0.5$ N are presented in Figures 14 and 15. The integral control
gain has be increased to $K_{I}=0.05$ to enhance PLL stabilization and
accelerate convergence in order to obtain a larger portion of the 1:3
resonance.

Figure 14: NFRC ( ) of the 1:3 resonance identified by PLLT ($F=0.5$N). The
starting point is the red marker. The reference solution is provided by HBM in
gray; light gray corresponds to unstable responses. The bifurcation points are
represented by circles. (a) Amplitude and (b) phase lag.

Figure 15: NFRC ( ) of the 1:2 resonance identified by PLLT ($F=0.5$N). The
starting point is the red marker. The reference solution is provided by HBM in
gray; light gray corresponds to unstable responses. The bifurcation points are
represented by circles. (a) Amplitude and (b) phase lag.

##  4 Experimental demonstration using a clamped-clamped beam

###  4.1 Experimental setup

The test sample is a thin straight aluminum beam of dimension
$2\mathrm{mm}\times 20\mathrm{mm}\times 1\mathrm{m}$, see Figure 16. Both ends
of the beam were clamped using two aluminum blocks bolted to a support frame
fixed on a vibration table. The clamped beam length could be adjusted by
changing the position of the aluminum blocks. A sliding groove of cross
section $0.2\mathrm{mm}\times 20\mathrm{mm}$ was made in the block to
prescribe the beam’s movements.

During the installation process, low initial tension or compression force
along the beam was ensured so as to avoid structural buckling. Such a slender
beam is expected to exhibit geometrically nonlinear behavior resulting in an
hardening effect due to mid-plane stretching. The first vibration mode is
considered in this study. The formulas given for secondary resonances in [6]
were useful for guiding experiment design.

An electromagnetic shaker (TIRA TV 51075) fed by a power amplifier (TIRA BAA
120) working in current mode was used to drive the beam into vibration. The
shaker was connected to the beam through a flexible stinger and an impedance
head (DYTRAN 5860B). The shaker was placed close to one clamped boundary in
order to minimize shaker–structure interaction. An accelerometer (DYTRAN
3035G) was placed closer to the beam center. The measured signals were
amplified by a signal conditioner (PCB PIEZOTRONICS 482C).

PLLT was realized using the real-time controller dSPACE MicroLabBox with the
control scheme implemented in Matlab/Simulink. The sampling frequency of the
fixed-step solver was set to be 10 kHz. An additional controller was used to
control the applied excitation toward a settled magnitude $F_{o}$ [21], as
schematized in Figure 17. The gains of the PI controller were both set to
0.01. The adaptive filter parameter $\mu$ was set to 0.001.

Despite its apparent simplicity, this set-up posed two important difficulties.
First, shaker-structure interaction remained an important issue for PPLT with
the consequence that it was not possible to find a single beam length for
which all the resonances of interest could be identified robustly. Second, the
set-up exhibited important day-to-day variability induced by the key role
played by room temperature on the dynamical properties of the beam.

Figure 16: Experimental setup.  Figure 17: Fundamental force controller.

###  4.2 Beam 1 (length=0.678m): primary and 1:2 subharmonic resonances

The NFRCs obtained from classical sine sweep testing (SST) and from PLLT
around the primary resonance as well as the backbone identified by PLLT with
$\Phi_{\mathrm{ref}}=\pi/2$ are displayed in Figure 18. Unlike the previous
section, the phase lag was measured between the acceleration and forcing
signals, so there is a difference of $\pi$ with respect to the results of
Section 3. The beam is seen to exhibit slight softening behavior at low
forcing levels followed by a more prominent hardening behavior at higher
levels. The softening behavior can be attributed to the clampings whereas the
hardening behavior is due to the geometrical nonlinearity. In addition to a
good correspondence between the backbone and the local maxima in the NFRCs
identified using PLLT, there is also good agreement between the NFRCs obtained
using PLLT and SST. The only difference is at $0.6$N where a jump is observed
for SST due to the folded nature of the NFRC. Conversely, PLLT is able to
unfold the NFRC and measure the complete branch including its unstable
portion. We note that, in Figure 18(b), the measured phase lag is not exactly
equal to $\pi/2$, which can be explained by the system’s variability between
different forcing levels as well as the existence of transients. However, the
relative error between $\hat{\Phi}_{1,1}$ and $\Phi_{\mathrm{ref}}$ remained
consistently below $2\%$.

Figure 18: NFRCs (solid lines in grayscale) and backbone (blue dashed line) of
the primary resonance identified by PLLT using $K_{P}=1$ and $K_{I}=2$ at
0.2N, 0.4N and 0.6N. The colored areas correspond to SST. (a) Amplitude and
(b) phase lag.

To verify the existence of the 1:2 subharmonic isola, an impact hammer was
used to trigger the jump from the main branch to the isola while the shaker
was exciting the structure at approximately twice the primary resonance
frequency. The isola was found to exist for $F\geq 2.3$N. The starting point
on the isola for forcing amplitudes between 2.5N and 2.7N was also obtained
using the impact hammer. The NFRCs measured with PLLT are depicted in Figure
19; the different harmonics measured at 2.7N are plotted in Figure 20. Large
PLL controller gains, i.e., $K_{P}=10$ and $K_{I}=25$, were necessary during
the identification process. From the starting point, the NFRCs were obtained
by sequentially decreasing the phase lag, resulting in an increase then
decrease in forcing frequency. This strategy allowed us to trace the high-
amplitude portion of the 1:2 isola around the turning point.

The backbone in Figure 19 was identified using $\Phi_{\mathrm{ref}}=0$, which
is not in accordance with the theoretical expectations for a Duffing-like type
of behavior. A possible explanation for this discrepancy is that the beam has
a non-negligible initial curvature resulting in a Helmholtz-like behavior for
which the resonant phase lag between the acceleration and forcing signals of a
1:2 subharmonic resonance is 0 or $\pi$ [36]. Nonetheless, the backbone is
found to provide a good approximation to the locus of the local maxima on the
1:2 isolas. Finally, the results of SST are also given in Figure 19. A
transition to the isola during SST was possible by applying perturbations. The
agreement between PLLT and SST is very good considering the challenges posed
by the inherent variability of the set-up and shaker-structure interactions.

Figure 19: NFRCs (solid lines in grayscale) and backbone (blue dashed line) of
the 1:2 subharmonic resonance identified by PLLT using $K_{P}=10$ and
$K_{I}=25$ at 2.5N, 2.6N and 2.7N. The colored areas correspond to SST. (a)
Amplitude and (b) phase lag. Figure 20: Different harmonics of the 1:2
subharmonic resonance measured at 2.7N: $\kappa=1$ (red), $\kappa=2$ (blue)
and $\kappa=3$ (yellow). The total amplitude is represented in black.

###  4.3 Beam 2 (length=0.754m): 2:1/3:1 superharmonic and 1:2 subharmonic
resonances

The results of the identification of the 3:1 superharmonic resonance are shown
in Figure 21. The backbone obtained by locking the phase at the theoretically-
expected value of $\pi/2$ for the acceleration ($-\pi/2$ for the displacement)
provides an excellent approximation of the local maxima of the three NFRCs.
The different harmonics measured at 4.5N in Figure 22 indicates that the third
harmonic in the vicinity of the frequency of the primary resonance is the
dominant harmonic.

Figure 21: NFRCs (solid lines in grayscale) and backbone (blue dashed line) of
the 3:1 superharmonic resonance identified by PLLT using $K_{P}=0.1$ and
$K_{I}=0.05$ at 3.5N, 4.0N and 4.5N. (a) Amplitude and (b) phase lag. Figure
22: Different harmonics of the 3:1 superharmonic resonance measured at 4.5N:
$\kappa=1$ (blue), $\kappa=3$ (red) and $\kappa=5$ (yellow). The total
amplitude is represented in black.

The identification of 2:1 superharmonic resonance was performed in a similar
manner. Despite that a variability of the resonance frequency was observed
during PLLT, a good correspondence between the NFRCs and the backbone curve
was obtained in Figure 23.

Figure 23: NFRCs (solid lines in grayscale) and backbone (blue dashed line) of
the 2:1 superharmonic resonance identified by PLLT using $K_{P}=1$ and
$K_{I}=2$ at 2.5N, 3.0N and 3.5N. (a) Amplitude and (b) phase lag.

Finally, the automatic transfer to the 1:2 isola and the subsequent
identification were also successfully achieved for beam 2. As shown in Figure
24, after some transient oscillations, the system was successfully attracted
by the 1:2 resonance and travelled along the isolated branch. Figure 24(d)
illustrates the growth of the half harmonic during the transfer from the main
branch to the subharmonic orbit.

Figure 24: Automatic transition to the detached 1:2 subharmonic resonance and
subsequent identification of the NFRCs and backbone. State transfer: $K_{P}=1$
and $K_{P}=0.5$; identification: $K_{P}=5$ and $K_{P}=20$.

##  5 Conclusions

In this study, a control-based vibration testing scheme for identifying both
primary and secondary resonances of nonlinear systems was proposed. The PLLT
scheme integrates digital adaptive filters into the phase detector module of
PLL. The experimental continuation of NFRCs exploited the fact that, unlike
the forcing frequency, the phase lag has a monotonic evolution around
nonlinear resonances. Thus, the NFRC can be obtained by sweeping the phase lag
at constant forcing amplitude. The experimental continuation of backbone
curves exploited the concept of a resonant phase lag. In this case, the
forcing amplitude is swept at constant phase lag.

The performance of PLLT was carefully assessed using numerical and physical
experiments. The NFRCs (with their stable and unstable portions) and backbones
of different superharmonic and subharmonic resonances were successfully
identified. A salient feature of PLLT is that both bifurcating and isolated
branches were reached from the rest position without any user interaction.
Throughout this study, the selection of the gains of the PI controller
remained a matter of trial and error. Further research concerning this aspect
is necessary.

## References

  * [1] G. Kerschen, K. Worden, A.F. Vakakis, and J.-C. Golinval.  Past, present and future of nonlinear system identification in structural dynamics.  Mechanical systems and signal processing, 20(3):505–592, 2006.
  * [2] D. Göge, M. Böswald, U. Füllekrug, and P. Lubrina.  Ground vibration testing of large aircraft–state-of-the-art and future perspectives.  In 25th International Modal Analysis Conference, 2007.
  * [3] J.-P. Noël, L. Renson, and G. Kerschen.  Complex dynamics of a nonlinear aerospace structure: experimental identification and modal interactions.  Journal of Sound and Vibration, 333(12):2588–2607, 2014.
  * [4] S. Schwarz, L. Kohlmann, A. Hartung, J. Gross, M. Scheel, and M. Krack.  Validation of a turbine blade component test with frictional contacts by phase-locked-loop and force-controlled measurements.  Journal of Engineering for Gas Turbines and Power, 142(5):051006, 2020.
  * [5] A. Opreni, A. Viazzacrro, A. Frangi, and C. Touze.  Model order reduction based on direct normal form: application to large finite element mems structures featuring internal resonance.  Nonlinear Dynamics, 105:1237–1272, 2021.
  * [6] A.H. Nayfeh and D.T. Mook.  Nonlinear oscillations.  John Wiley & Sons, 2008.
  * [7] R.M. Rosenberg.  Normal modes of nonlinear dual-mode systems.  Journal of Applied Mechanics, 27(2):263–268, 1960.
  * [8] R.M. Rosenberg.  The normal modes of nonlinear n-degree-of-freedom systems.  Journal of Applied Mechanics, 29(1):7, 1962.
  * [9] G. Kerschen, M. Peeters, J.-C. Golinval, and A.F. Vakakis.  Nonlinear normal modes, part i: A useful framework for the structural dynamicist.  Mechanical systems and signal processing, 23(1):170–194, 2009.
  * [10] S. Shaw and C. Pierre.  Non-linear normal modes and invariant manifolds.  Journal of sound and Vibration, 150(1):170–173, 1991.
  * [11] S. Shaw and C. Pierre.  Normal modes for non-linear vibratory systems.  Journal of sound and vibration, 164(1):85–124, 1993.
  * [12] G. Haller and S. Ponsioen.  Nonlinear normal modes and spectral submanifolds: Existence, uniqueness and use in model reduction.  Nonblinear Dynamics, 86:1493–1534, 2016.
  * [13] T.L. Hill, A. Cammarano, S.A. Neild, and D.J. Wagg.  Interpreting the forced responses of a two-degree-of-freedom nonlinear oscillator using backbone curves.  Journal of Sound and Vibration, 349:276–288, 2015.
  * [14] C. Touzé, A. Vizzaccaro, and O. Thomas.  Model order reduction methods for geometrically nonlinear structures: a review of nonlinear techniques.  Nonlinear Dynamics, 105(2):1141–1190, 2021.
  * [15] M. Peeters, G. Kerschen, and J.-C. Golinval.  Dynamic testing of nonlinear vibrating structures using nonlinear normal modes.  Journal of Sound and Vibration, 330(3):486–509, 2011.
  * [16] M. Peeters, G. Kerschen, and J.-C. Golinval.  Modal testing of nonlinear vibrating structures based on nonlinear normal modes: Experimental demonstration.  Mechanical Systems and Signal Processing, 25(4):1227–1247, 2011.
  * [17] S. Peter and R.I. Leine.  Excitation power quantities in phase resonance testing of nonlinear systems with phase-locked-loop excitation.  Mechanical Systems and Signal Processing, 96:139–158, 2017.
  * [18] M. Scheel, S. Peter, R.I. Leine, and M. Krack.  A phase resonance approach for modal testing of structures with nonlinear dissipation.  Journal of Sound and Vibration, 435:56–73, 2018.
  * [19] V. Denis, M. Jossic, C. Giraud-Audine, B. Chomette, A. Renault, and O. Thomas.  Identification of nonlinear modes using phase-locked-loop experimental continuation and normal form.  Mechanical Systems and Signal Processing, 106:430–452, 2018.
  * [20] A. Givois, J.-J. Tan, C. Touzé, and O. Thomas.  Backbone curves of coupled cubic oscillators in one-to-one internal resonance: bifurcation scenario, measurements and parameter identification.  Meccanica, 55(3):481–503, 2020.
  * [21] G. Abeloos, F. Müller, E. Ferhatoglu, M. Scheel, C. Collette, G. Kerschen, M.R.W. Brake, P. Tiso, L. Renson, and M. Krack.  A consistency analysis of phase-locked-loop testing and control-based continuation for a geometrically nonlinear frictional system.  Mechanical Systems and Signal Processing, 170:108820, 2022.
  * [22] U. Parlitz and W. Lauterborn.  Superstructure in the bifurcation set of the duffing equation.  Physics Letter A, 107(8):351–355, 1985.
  * [23] A. Vizzaccaro, G. Gobat, A. Frangi, and C. Touzé.  Direct parametrisation of invariant manifolds for generic non-autonomous systems including superharmonic resonances.  arXiv preprint arXiv:2306.09860, 2023.
  * [24] A.P. Chauhan and D.G. Ashwell.  On the low-frequency drumming of slightly bowed structures.  Journal of Sound and Vibration, 14(4):475–489, 1971.
  * [25] N. Yamaki, K. Otomo, and M. Chiba.  Non-linear vibrations of a clamped circular plate with initial deflection and initial edge displacement, part ii: Experiment.  Journal of Sound Vibration, 79(1):43–59, 1981.
  * [26] N. Yamaki, K. Otomo, and M. Chiba.  Nonlinear vibrations of a clamped rectangular plate with initial deflection and initial edge displacement—part ii: Experiment.  Thin-Walled Structures, 1(2):101–119, 1983.
  * [27] E.A. Neppiras.  Subharmonic and other low-frequency emission from bubbles in sound-irradiated liquids.  The Journal of the Acoustical Society of America, 46(3B):587–601, 1969.
  * [28] H.J. Janssen, L. Beerden, and E.L.M. Flerackers.  An experimental look at the resonant behaviour of a nonlinear lc circuit.  European Journal of Physics, 5(2):94, 1984.
  * [29] H. Lin, S.C.S. Yim, and O. Gottlieb.  Experimental investigation of response stability and transition behaviour of a nonlinear ocean structural system.  Ocean engineering, 25(4-5):323–343, 1998.
  * [30] H.S. Husain and F. Hussain.  Experiments on subharmonic resonance in a shear layer.  Journal of Fluid mechanics, 304:343–372, 1995.
  * [31] M. Zhang, Y. Shen, L. Xiao, and W. Qu.  Application of subharmonic resonance for the detection of bolted joint looseness.  Nonlinear dynamics, 88:1643–1653, 2017.
  * [32] M. Volvert and G. Kerschen.  Phase resonance nonlinear modes of mechanical systems.  Journal of Sound and Vibration, 511:116355, 2021.
  * [33] G. Abeloos, L. Renson, C. Collette, and G. Kerschen.  Stepped and swept control-based continuation using adaptive filtering.  Nonlinear Dynamics, 104(4):3793–3808, 2021.
  * [34] G. Abeloos.  Control-based methods for the identification of nonlinear structures.  PhD thesis, University of Liège, 2022.
  * [35] S. Haykin.  Adaptive filter theory.  Prentice Hall, 1996.
  * [36] M. Volvert.  Resonant phase lags of nonlinear mechanical systems.  PhD thesis, University of Liège, 2024.
  * [37] T. Detroux, L. Renson, L. Masset, and G. Kerschen.  The harmonic balance method for bifurcation analysis of large-scale nonlinear mechanical systems.  Computer Methods in Applied Mechanics and Engineering, 296:18–38, 2015.