# Leveraging Geometrical Acoustic Simulations of Spatial Room Impulse
Responses for Improved Sound Event Detection and Localization

###### Abstract

As deeper and more complex models are developed for the task of sound event
localization and detection (SELD), the demand for annotated spatial audio data
continues to increase. Annotating field recordings with 360∘ video takes many
hours from trained annotators, while recording events within motion-tracked
laboratories are bounded by cost and expertise. Because of this, localization
models rely on a relatively limited amount of spatial audio data in the form
of spatial room impulse response (SRIR) datasets, which limits the progress of
increasingly deep neural network based approaches. In this work, we
demonstrate that simulated geometrical acoustics can provide an appealing
solution to this problem. We use simulated geometrical acoustics to generate a
novel SRIR dataset that can train a SELD model to provide similar performance
to that of a real SRIR dataset. Furthermore, we demonstrate using simulated
data to augment existing datasets, improving on benchmarks set by state of the
art SELD models. We explore the potential and limitations of geometric
acoustic simulation for localization and event detection. We also propose
further studies to verify the limitations of this method, as well as further
methods to generate synthetic data for SELD tasks without the need to record
more data.

Index Terms—  Acoustic Simulation, Localization, Data Augmentation

##  1 Introduction

Sound event detection and localization (SELD) is the union of two active
fields of research; sound event detection (SED), and localization, or
direction-of-arrival (DoA) estimation. Expanding the scene-description
capabilities of SED with the spatiotemporal characterization of localization
sees applications ranging from autonomous robot navigation [1] and urban
monitoring [2], to speaker diarization [3] and immersive experiences in
virtual and augmented reality devices.

While earlier techniques for SELD have focused on traditional signal-
processing or parametric models such as [4, 5, 6, 7], recent literature is
dominated by deep neural network (DNN) approaches, which have been shown
superior performance in both pure localization [8, 9] as well as joint SELD
tasks [10]. A surge of interest in this field can be attributed to the
introduction of SELD as a task in the DCASE2019 challenge [11]. This challenge
included the release of several 4-channel audio datasets with spatial and
temporal annotations for sound events. The audio was generated by convolving
sound events with spatial room impulse responses (SRIRs) recorded in 5
separate rooms at 504 unique azimuth-elevation-distance combinations. This was
further iterated upon by the SELD challenge in DCASE2020 [12] with 13 unique
rooms. Recently, DCASE2022 reintroduced the challenge with hand-annotated
real-world recordings in the STARSS22 dataset [13], providing one of the first
datasets with real-world data upon which to evaluate SELD models. This dataset
uses a combination of 360∘ video, and motion capture to extract spatiotemporal
annotations that were manually validated. In addition to this, the DCASE2022
dataset also included a release of the SRIRs used to generate the training
data for the SELD task, as well as the code for the generator itself, allowing
users to generate their own annotated spatial data [14]. This data included
SRIRs measured over a wide range of positions over 9 different rooms in
Tampere University’s campus.

These datasets are unique in the density of SRIR measurements across
particular paths, the variety of acoustic enclosures, and the large amount of
SRIRs. Because of their scale, visibility, and quality, these datasets have
become some of the most cited SRIR datasets for DNN-based approaches to SELD,
because of their ability to meet the data requirements of these highly-
parametrized models. Despite this, these datasets are still severely limited
by the recording procedure of SRIRs, which require time, expertise, and a low-
noise environment to produce at a high quality. Increasing the spatial
density, variety of trajectories types, and number of trajectory paths becomes
multiplicatively time consuming to develop. Furthermore, the range of rooms in
which these measurements can be recorded is inherently limited by the
limitations of the recording facilities, usually limited to a dozen rooms or
so in the best of cases. However, without a wide range of acoustic
environments to perform these measurements, generalization to a variety of
unseen acoustic environments becomes impossible.

Physical acoustic simulations provide an attractive solution to the
limitations of field-recorded SRIRs. Acoustic simulation is typically split
into two categories: wave based methods, which simulate the propagation of
sound waves through physical media, and geometric modeling methods, which
model the transportation of acoustic energy through acoustic rays, mimicking
popular methods for modeling optical rays. Geometrical acoustics approximate
the wavelength of the propagating sound to have wavelength relatively small
compared to the room geometries of interest, and neglects wave effects such as
diffraction or scattering. Nevertheless because of ease of implementation and
computational efficiency, geometrical acoustic modeling methods have seen wide
success in several tasks, including modeling architectural acoustics [15] and
room parameter estimation [16].

In this work, we propose utilizing one method of geometrical acoustics
modeling, the image source method, to generate simulated SRIRs for training
DNN models for SELD. We demonstrate the effectiveness of this simulated method
for SRIR generation using the framework and data provided in previous DCASE
SELD challenges. By creating an audio dataset from simulated SRIRs, we train a
SELD model with similar performance to one utilizing real-world SRIRs. By
directly compares simulated SRIRs with a datasets of recorded SRIRs of the
similar size, room geometries, and DoA distributions, we demostrate the
downstream effects of simulation in place of recording as being relatively
minimal, differing our work from prior studies [17]. Furthermore, we augment a
typical SRIR dataset with simulated SRIRs, training models that outperform
those trained solely on recorded SRIRs. Finally we propose further experiments
to explore the use of simulated SRIRs for training SELD models. The code
associated with this work is released in an open-source github repository
111https://github.com/ChrisIck/DCASE_Synth_Data to further work in using
synthetic SRIRs for training DNN models.

##  2 Acoustic Simulations

###  2.1 The Image Source Method

The image source method (ISM) is a technique used in architectural acoustics
and room modeling to predict the sound field in enclosed spaces [18]. The ISM
considers the primary sound source and virtual images reflected by the room’s
boundaries. These virtual sources are assumed to emit sound with the same
magnitude and phase as the primary source, but with a delay due to the
additional path length traveled. Typically, this starts with defining the room
geometry, including the positions and shapes of the walls, ceiling, and floor.
For each reflecting surface, virtual image sources are “mirrored” across the
boundary. The number of virtual sources depends on the order of reflections
considered. From here, the interaction between the primary sound source and
the virtual image sources is calculated by determining the path lengths, time
delays, and attenuation factors associated with each source-receiver
combination, as well as the material properties of each surface through which
the sound path is reflected upon. By summing the contributions of the primary
source and its image sources, the sound field at various locations in the room
can be predicted, providing an estimation of the sound pressure level, arrival
times, and directivity patterns.

It’s important to note several limitations of the ISM model. The ISM
implicitly assumes that all surfaces are perfectly reflective and flat, with
idealized acoustic properties, which fails to account for acoustic effects
such as scattering or diffraction. Furthermore, the order of
reflections/virtual sources scales the computational cost exponentially,
meaning compute late-stage reflections in an SRIR prohibitively expensive.

Despite its limitations, the image source method is widely used due to its
computational efficiency and effectiveness in predicting sound fields in
enclosed spaces. Because localization is more reliant on direct sounds/early
reflections, the limitations caused by use of the the ISM for computing SRIRs
can be expected to be relatively minimal.

###  2.2 The TAU-SRIR Dataset

To validate the use of ISM-generated SRIRs in a direct-comparison, we take the
existing TAU-SRIR database [14] as an example database for which well
established metrics for SELD have been measured.

The TAU-SRIR database contains SRIRs recorded in 9 different rooms throughout
Tampere University’s campus. Each SRIR is computed by recording a maximum-
length sequence (MLS) played through a loudspeaker, recorded on an Eigenmike
spherical microphone array. Each SRIR was downsampled to 24kHz and truncated
at 300ms, resulting in $7200$ samples per RIR. The data is stored in a
4-channel audio corresponding to a tetrahedral microphone array with the
geometry in spherical coordinates $(\phi,\theta)$, specified in Table 1. For
each room, the position of the microphone array was provided.

|  Azimuth ($\phi$) |  Elevation ($\theta$)
---|---|---
M1 | $45^{\circ}$ | $35^{\circ}$
M2 | $-45^{\circ}$ | $-35^{\circ}$
M3 | $135^{\circ}$ | $-35^{\circ}$
M4 | $-135^{\circ}$ | $35^{\circ}$
Table 1: Microphone Geometry for TAU-SRIR dataset. Each microphone is 4.2cm
from the center, and is modeled with a hypercardioid response.

SRIRs were measured along either circular or linear traces at fixed distance
from the microphone array along the z-axis at a number of trajectory groups,
separated by distance and reflection across the axis of the microphone array
in the case of linear traces. Circular trajectory groups had a specified
radius of orbit, whereas linear trajectories had a specified start and end
point in 3D space. Each trajectory was repeated at a number of different
heights, and each trajectory had a fixed number SRIR measurements and
corresponding DoA measurements recorded as Cartesian components of a unit
vector. The number of SRIR measurements vary across different
trajectories/heights, spaced in roughly $1^{\circ}$ increments. The total
number measurements can be seen in Table 2

Room Name | Traj. type | $N_{t}$ | $N_{h}$ | $N_{\text{SRIRs}}$
---|---|---|---|---
Bomb shelter | Circular | 2 | 9 | 6480
Gym | Circular | 2 | 9 | 6480
PB132 | Circular | 2 | 9 | 6480
PC226 | Circular | 2 | 9 | 6480
SA203 | Linear | 6 | 3 | 1594
SC203 | Linear | 4 | 5 | 1592
SE203 | Linear | 4 | 4 | 1760
TB103 | Linear | 4 | 3 | 1184
TC352 | Circular | 2 | 9 | 6480
Table 2: Trajectory information for rooms contained in the TAU-SRIR dataset
[14]. Each room contains trajectories across a number of trajectory groups
($N_{t}$) and a number of heights ($N_{h}$), for a total of $N_{t}\times
N_{h}$ trajectories per room. Each trajectory is sampled in roughly
$1^{\circ}$ increments.

###  2.3 Room Simulation

Figure 1: The log-scale energy decay from a sample of RIRs from a singular
room. The region in the dashed lines is roughly linear, suggesting it
corresponds to mid-late reflections, and is used for an estimate of RT60

We recreate this dataset using the python package pyroomacoustics [19], a
pythonic implementation of the ISM, that has demonstrated use in
implementations of various algorithms for beamforming, direction finding,
adaptive filtering, source separation, and single channel denoising.

To replicate the acoustic conditions of each of the rooms in the TAU-SRIR
dataset, we randomly sampled the RIRs uniformly in each room until we had a
sample of 5 single-channel RIRs. Using the Schroeder method [20], we estimated
the RT60 of each room by hand-selecting the early decay of the energy-decay
function of the RIR samples and computing the linear fit. This can be see in
in Figure 1. Using the inverse Sabine formula, we used this to estimate the
mean absorption coefficient of the rooms and the number of required reflection
orders to approximate a room of a similar RT60. We combined these parameters
with the geometry estimations from the TAU-SRIR dataset to construct virtual
rooms matching those of the 9 rooms in the TAU-SRIR dataset. To this room, we
added a virtual tetrahedral microphone with the geometry specified in Table 1,
with each virtual microphone using a hypercardioid response pattern centered
at the position specified in the TAU-SRIR dataset.

Figure 2: SRIR measurement positions reconstructed from the DoAs provided in
the TAU-SRIR dataset (blue), compared to the path specified by the height,
radius, and position labels provided (red)

To estimate the positions of the SRIRs in space corresponding to the DoA
measurements in the TAU-SRIR dataset, we chose points along the DoA that most
closely matched the corresponding path on the line specified by the trajectory
in the dataset. Circular datasets specified their height and radius, and were
centered along the same z-line as the microphone array. Linear trajectories
specified the start point and end points of their traces. These matched points
were estimated by projecting the DoA vectors onto a cylinder that matched the
radius of the trajectory, after translation by the position of the microphone
array and the height of the trajectory of interest (see Figure 2). Once these
points were estimated, they were placed into the simulated room in a position
relative to the virtual microphone array, and the 4-channel SRIR was computed
with the ISM for each point, at the sample rate of 2400kHz. To match the
dimensions of the original TAU-SRIR dataset, we truncate the SRIRs to 300ms,
providing us with an SRIR with dimensions of $7200\times 4$ for each point.

##  3 Methodology

To evaluate the performance of this dataset in SELD tasks, we generated a
dataset of audio events consistent with the methodology of dataset generation
for training the baseline SELD model in the DCASE 2022 challenge [21]. We
generated 3 datasets, one using the original SRIR database, which we will
refer to as the TAU-SRIR datset. We also generated a dataset using only the
synthetic SRIRs, which we refer to as the SIM-SRIR dataset. Finally, we
generated a third dataset that equally samples both the original and simulated
SRIRs, which we will refer to as the augmented SRIR dataset, or AUG-SRIR.

###  3.1 Data Generation

To generate our annotated spatialized audio, we followed the procedure used in
DCASE2019-2021, by convolving various sound events with SRIRs.

The sound events were drawn from the FSD50k audioset [22], a subset containing
over 20k sound events of 13 classes selected for the DCASE challenge. These
sound events were spatialized into virtual recordings, each corresponds to a
singular room, allowing for up to three concurrently active sources. The
sources can be static or dynamic, with equal probability, and the dynamic
sources can move at slow ($10^{\circ}$/sec), moderate ($20^{\circ}$/sec), or
fast ($40^{\circ}$/sec) angular speeds. Each sample lasted 60 seconds, 40 of
which had at least one active event class.

For each of the 3 SRIR datasets, 1200 recordings were created in separate
folds, 900 for training and 300 for validation. The training and validation
sets used 6 and 3 rooms respectively, such that none of the same rooms
overlapped both folds.

###  3.2 Model

The model architecture is identical to the one used in the DCASE2022 Task 3
challenge baseline [13]; a SELDnet style CRNN with multi-ACCDOA representation
for co-occuring events [23]. The model takes $T$ frames of an STFT time-
frequency representation of the multichannel features, and outputs $T/5\times
N\times C\times 3$ vector coordinates, and $N$ is the assumed maximum co-
occurring events, in our case $3$.

The input features are 4-channel 64-band log-mel spectrograms combined with
SALSA-lite spatial features, all of which are truncated to include bins up to
9kHz, without mel-band aggregation following [24]).

Each model was trained on the training folds generated from the SRIR datasets
described above. In addition to this, data from the Sony-TAu Realistic Spatial
Soundscapes 2022 (STARSS22) [13] was added for training, using the 54
development sound mixtures for training, but witholding the remaining 52 clips
for evaluation of the models, ensuring the results were exclusively on real
recorded data from unseen rooms. The models were trained for 100 epochs each,

###  3.3 Evaluation

Evaluation was completed using join localization-detection metrics established
in the DCASE 2020 challenge. The detection metrics used were error-rate and F1
score for a spatial threshold within $20^{\circ}$ ($ER_{20^{\circ}}$ and
$F_{20^{\circ}}$). F1 score was macro-averaged to account for class
distribution differences in the FSD50k audio subset used. Localization metrics
are class dependent localization error and recall ($LE_{\text{CD}}$ and
$LR_{\text{CD}}$).

##  4 Results

Despite the coarse physical approximations made by the ISM, the entirely
synthetic SRIRs generated with this process performed nearly as well as the
SRIRs recorded in real world settings. Furthermore, the dataset of real SRIRs
augmented by synthetic SRIRs outperformed both by a narrow margin, showing
benefits of geometrical acoustics simulation for data augmentation for SELD
tasks.

Regarding the cross-class average performance of the models in Table 3, we can
see that for our classification metrics, all three models perform relatively
similarly, with a slight performance edge to the AUG-SRIR dataset trained
models. Looking into the per-class results in figure 3, we can see that
generally, all three models struggle with similar classes (telephone,
laughter, door), but the AUG-SRIR dataset outperforms both in certain classes
for which both other models perform poorly on (Water tap/Faucet, and Knock).

Looking at the localization based results, it appears that some amount of the
performance differences between the SIM-SRIR trained models and the TAU-SRIR
trained models can possibly be attributed to model fine-tuning; while SIM-SRIR
trained models had poor localization performance on certain results (Water
tap/Faucet, and Knock), the AUG-SRIR model outperformed the baseline TAU-SRIR
dataset. This suggests that the SIM-SRIR datasets are actually providing
beneficial information for these sound classes missing from the TAU-SRIR
datasets. With more thorough model tuning, it’s possible that the performance
for SIM-SRIR trained models it even closer to that of the baseline.

| $ER_{20^{\circ}}$ | $F_{20^{\circ}}$ | $LE_{\text{CD}}$ | $LR_{\text{CD}}$
---|---|---|---|---
TAU-SRIR | 0.71 | 14.4% |  55.1∘ | 39.2%
SIM-SRIR | 0.73 | 13.0% |  79.6∘ | 34.8%
AUG-SRIR | 0.75 |  16.3% | 52.3∘ | 42.3%
Table 3: The cross-class evaluation metrics for models trained on data
generated from different SRIR datasets Figure 3: Per-class results of models
trained on each dataset for F-measure within $20^{\circ}$, localization error,
and localization recall.

##  5 Conclusion

In this work we demonstrated the potential of using acoustic simulation to
generate spatial audio data for training SELD models. We’ve shown that
simulated SRIR data can improve the performance of SELD models as a form of
data augmentation. In addition to this, we’ve shown that simulated SRIRs,
while not as effective as those recorded in real acoustic environments, can be
used to effectively train SELD models, removing the relatively high cost of
producing additional data for similarly performing results in a relatively
limited setting. Generating larger volumes of SRIRs over a wider range of
acoustic conditions could provide even better results than these baselines,
potentially demonstrating greater robustness over varying acoustic
environments. Furthermore, using a high-volume of simulated SRIRs to train a
model, and using a hold-out of limited high-quality real-world data to fine
the model could produce SoTA results. This result is promising for future
experiments involve SRIRs for use in acoustic simulation data. Understanding
the requirements for angular density in dynamic SRIR recordings can help
inform future dataset collection practices, as well as the robustness of these
models to noise; limited work was done exploring the effect of noise on the
models trained with simulated SRIRs. Further ablation studies are necessary to
understand the limitations of geometrical acoustic methods for SELD-based
tasks, but these early experiments suggest that these can provide a low-
resource alternative to real-world SRIR recordings.

## References

  * [1] L. Bondi, G. Chuang, C. Ick, A. Dave, C. Shelton, B. Coltin, T. Smith, and S. Das, “Acoustic imaging aboard the international space station (iss): Challenges and preliminary results,” in _ICASSP 2022 - 2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2022, pp. 5108–5112.
  * [2] F. Stevens and D. Murphy, “Spatial impulse response measurement in an urban environment,” in _Audio Engineering Society Conference: 55th International Conference: Spatial Audio_ , Aug 2014. [Online]. Available: http://www.aes.org/e-lib/browse.cfm?elib=17355
  * [3] J. Wang, Y. Liu, B. Wang, Y. Zhi, S. Li, S. Xia, J. Zhang, F. Tong, L. Li, and Q. Hong, “Spatial-aware speaker diarizaiton for multi-channel multi-party meeting,” in _Interspeech 2022_. ISCA, sep 2022. [Online]. Available: https://doi.org/10.21437%2Finterspeech.2022-11412
  * [4] T. Butko, F. G. Pla, C. Segura, C. Nadeu, and J. Hernando, “Two-source acoustic event detection and localization: Online implementation in a smart-room,” in _2011 19th European Signal Processing Conference_ , 2011, pp. 1317–1321.
  * [5] C. J. Grobler, C. P. Kruger, B. J. Silva, and G. P. Hancke, “Sound based localization and identification in industrial environments,” in _IECON 2017 - 43rd Annual Conference of the IEEE Industrial Electronics Society_ , 2017, pp. 6119–6124.
  * [6] R. Chakraborty and C. Nadeu, “Sound-model-based acoustic source localization using distributed microphone arrays,” in _2014 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2014, pp. 619–623.
  * [7] K. Lopatka, J. Kotus, and A. Czyzewski, “Detection, classification and localization of acoustic events in the presence of background noise for acoustic surveillance of hazardous situations,” _Multimedia Tools Appl._ , vol. 75, no. 17, p. 10407–10439, sep 2016. [Online]. Available: https://doi.org/10.1007/s11042-015-3105-4
  * [8] S. Adavanne, A. Politis, and T. Virtanen, “Direction of arrival estimation for multiple sound sources using convolutional recurrent neural network,” _CoRR_ , vol. abs/1710.10059, 2017. [Online]. Available: http://arxiv.org/abs/1710.10059
  * [9] L. Perotin, R. Serizel, E. Vincent, and A. Guérin, “Crnn-based multiple doa estimation using acoustic intensity features for ambisonics recordings,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 13, no. 1, pp. 22–33, 2019.
  * [10] S. Adavanne, A. Politis, J. Nikunen, and T. Virtanen, “Sound event localization and detection of overlapping sources using convolutional recurrent neural networks,” _CoRR_ , vol. abs/1807.00129, 2018. [Online]. Available: http://arxiv.org/abs/1807.00129
  * [11] S. Adavanne, A. Politis, and T. Virtanen, “Localization, detection and tracking of multiple moving sound sources with a convolutional recurrent neural network,” _CoRR_ , vol. abs/1904.12769, 2019. [Online]. Available: http://arxiv.org/abs/1904.12769
  * [12] A. Politis, S. Adavanne, and T. Virtanen, “A dataset of reverberant spatial sound scenes with moving sources for sound event localization and detection,” in _Proceedings of the Detection and Classification of Acoustic Scenes and Events 2020 Workshop (DCASE2020)_ , Tokyo, Japan, November 2020, pp. 165–169. [Online]. Available: https://dcase.community/workshop2020/proceedings
  * [13] A. Politis, K. Shimada, P. Sudarsanam, S. Adavanne, D. Krause, Y. Koyama, N. Takahashi, S. Takahashi, Y. Mitsufuji, and T. Virtanen, “STARSS22: A dataset of spatial recordings of real scenes with spatiotemporal annotations of sound events,” in _Proceedings of the 8th Detection and Classification of Acoustic Scenes and Events 2022 Workshop (DCASE2022)_ , Nancy, France, November 2022, pp. 125–129. [Online]. Available: https://dcase.community/workshop2022/proceedings
  * [14] A. Politis, S. Adavanne, and T. Virtanen, “TAU Spatial Room Impulse Response Database (TAU- SRIR DB),” Apr. 2022. [Online]. Available: https://doi.org/10.5281/zenodo.6408611
  * [15] J. B. Allen and D. A. Berkley, “Image method for efficiently simulating small‐room acoustics,” _The Journal of the Acoustical Society of America_ , vol. 65, no. 4, pp. 943–950, 04 1979. [Online]. Available: https://doi.org/10.1121/1.382599
  * [16] C. Ick, A. Mehrabi, and W. Jin, “Blind acoustic room parameter estimation using phase features,” in _ICASSP 2023 - 2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2023, pp. 1–5.
  * [17] F. B. Gelderblom, Y. Liu, J. Kvam, and T. A. Myrvoll, “Synthetic data for dnn-based doa estimation of indoor speech,” in _ICASSP 2021 - 2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2021, pp. 4390–4394.
  * [18] J. B. Allen and D. Berkley, “Image method for efficiently simulating small‐room acoustics,” _The Journal of the Acoustical Society of America_ , vol. 65, no. 4, pp. 943–950, 04 1979. [Online]. Available: https://doi.org/10.1121/1.382599
  * [19] R. Scheibler, E. Bezzam, and I. Dokmanic, “Pyroomacoustics: A python package for audio room simulations and array processing algorithms,” _CoRR_ , vol. abs/1710.04196, 2017. [Online]. Available: http://arxiv.org/abs/1710.04196
  * [20] M. R. Schroeder, “New method of measuring reverberation time,” _The Journal of the Acoustical Society of America_ , vol. 37, no. 3, pp. 409–412, 1965\. [Online]. Available: https://doi.org/10.1121/1.1909343
  * [21] K. Shimada, Y. Koyama, S. Takahashi, N. Takahashi, E. Tsunoo, and Y. Mitsufuji, “Multi-accdoa: Localizing and detecting overlapping sounds from the same class with auxiliary duplicating permutation invariant training,” in _IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , Singapore, Singapore, May 2022.
  * [22] E. Fonseca, X. Favory, J. Pons, F. Font, and X. Serra, “Fsd50k,” Oct. 2020. [Online]. Available: https://doi.org/10.5281/zenodo.4060432
  * [23] K. Shimada, Y. Koyama, S. Takahashi, N. Takahashi, E. Tsunoo, and Y. Mitsufuji, “Multi-accdoa: Localizing and detecting overlapping sounds from the same class with auxiliary duplicating permutation invariant training,” 2022.
  * [24] T. N. T. Nguyen, D. L. Jones, K. N. Watcharasupat, H. Phan, and W.-S. Gan, “SALSA-lite: A fast and effective feature for polyphonic sound event localization and detection with microphone arrays,” in _ICASSP 2022 - 2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, may 2022\. [Online]. Available: https://doi.org/10.1109%2Ficassp43922.2022.9746132


# Self-Excited Dynamics of Discrete-Time Lur’e Systems

Juan A. Paredes<EMAIL_ADDRESS>Syed Aseem Ul Islam<EMAIL_ADDRESS>Omran Kouba<EMAIL_ADDRESS>Dennis S. Bernstein<EMAIL_ADDRESS>University of Michigan, Ann Arbor, Michigan, 48109, USA  Higher Institute for
Applied Sciences and Technology, Damascus, Syria

###### Abstract

Self-excited systems arise in numerous applications, such as biochemical
systems, fluid-structure interaction, and combustion. This paper analyzes a
discrete-time Lur’e system with a piecewise-linear saturation feedback
nonlinearity. The main result provides sufficient conditions under which the
Lur’e system is self-excited in the sense that its response is bounded and
nonconvergent.

###### keywords:

self-oscillation , self-excitation , discrete-time , nonlinear feedback ,
Lur’e system

##  1 Introduction

A self-excited system has the property that the input is constant but the
response is oscillatory. Self-excited systems arise in numerous applications,
such as biochemical systems, fluid-structure interaction, and combustion. The
classical example of a self-excited system is the van der Pol oscillator,
which has two states whose asymptotic response converges to a limit cycle. A
self-excited system, however, may have an arbitrary number of states and need
not possess a limit cycle. Overviews of self-excited systems are given in [1,
2].

Models of self-excited systems are typically based on the relevant physics of
the application. From a systems perspective, the main interest is in
understanding the features of the components of the system that give rise to
self-sustained oscillations. Understanding these mechanisms can illuminate the
relevant physics in specific domains and provide unity across different
domains.

A unifying model for self-excited systems is a feedback loop involving linear
and nonlinear elements; systems of this type are called Lur’e systems. Lur’e
systems have been widely studied in the classical literature on stability
theory [3]. Within the context of self-excited systems, Lur’e systems under
various assumptions are considered in [4, 5, 6, 7, 8, 9, 1, 10, 11, 12]. Self-
oscillating discrete-time systems are considered in [13, 14].

Roughly speaking, self-excited oscillations arise from a combination of
stabilizing and destabilizing effects. Destabilization at small signal levels
causes the response to grow from the vicinity of an equilibrium, whereas
stabilization at large signal levels causes the response to decay when the
state is far from an equilibrium. In particular, negative damping at low
signal levels and positive damping at high signal levels is the mechanism that
gives rise to a limit cycle in the van der Pol oscillator [15, pp. 103–107].
Note that, although systems with limit-cycle oscillations are self-excited,
the converse need not be true since the response of a self-excited system may
oscillate without the trajectory reaching a limit cycle. Alternative
mechanisms exist, however; for example, time delays are destabilizing, and
Lur’e models with time delay have been considered as models of self-excited
systems [16].

The present paper considers a discrete-time Lur’e system with asymptotically
stable linear dynamics, a zero at 1, and a piecewise-linear saturation
feedback nonlinearity. For this Lur’e system, sufficiently large scalings of
the asymptotically stable dynamics yield closed-loop unstable dynamics while
the saturation function operates in its linear region. Under large signal
levels, the saturation function yields a constant signal, which effectively
breaks the loop, thus allowing the open-loop dynamics to stabilize the
response. The zero at 1 acts as a high-pass filter, which ensures that the
response does not converge, whereas the saturation function yields a constant
signal. Hence, while the saturation function operates in the nonsaturated
region, the closed-loop system is unstable, and, while the saturation function
operates in the saturated region and yields a constant signal, the closed-loop
system is asymptotically stable and has a constant input. The contribution of
the present paper is to prove that this model structure yields self-excited
oscillations for sufficiently large scalings of the asymptotically stable
dynamics. A preliminary study of self-excited oscillations in a similar
discrete-time Lur’e model was performed in [17]. However, the present paper
goes far beyond [17] in breadth and depth of the analysis of these systems.

The contents of the paper are as follows. Section 2 considers a simple
discrete-time linear feedback system and analyzes its transfer function to
study the range of values of the linear dynamics scalings for which the
closed-loop system is asymptotically stable. Section 3 considers the same
linear feedback system and analyzes its state space model to study the
conditions under which the response of the closed-loop system does not
converge and is not bounded. Section 4 extends the problem in Sections 2 and 3
by including a saturation nonlinearity. Under certain conditions, this
discrete-time Lur’e system is shown to have a bounded, non-convergent response
for sufficiently large values of the loop gain. Section 5 presents an example
that illustrates the conditions for self-excitation presented in Section 4.

Nomenclature.
${\mathbb{R}}\stackrel{{\scriptstyle\triangle}}{{=}}(-\infty,\infty),$
${\mathbb{C}}$ denotes the complex numbers, ${\mathcal{R}}$ denotes range,
${\mathcal{N}}$ denotes null space, $\overline{(\cdot)}$ denotes complex
conjugate, $(\cdot)^{*}$ denotes complex conjugate transpose, $\|\cdot\|$
denotes the Euclidean norm on ${\mathbb{C}}^{n},$ and ${\bf q}$ denotes the
forward shift operator. For a polynomial $p,$ define
$\operatorname{spr}(p)\stackrel{{\scriptstyle\triangle}}{{=}}\max\\{|z|\colon
z\in{\mathbb{C}}\mbox{ and }p(z)=0\\}$, and for a matrix
$A\in{\mathbb{R}}^{n\times n},$ let $\operatorname{spec}(A)$ denote the set of
eigenvalues of $A,$ let $\chi_{A}$ denote the characteristic polynomial of
$A,$ and let $\operatorname{spr}(A)$ denote the spectral radius of $A.$ For
$\gamma>0,$ ${\rm sat}_{\gamma}$ denotes the saturation function, such that,
for all $x\in[-\gamma,\gamma],$ ${\rm sat}_{\gamma}(x)=x,$ and, for all
$|x|>\gamma,$ ${\rm sat}_{\gamma}(x)=(\operatorname{sign}x)\gamma.$

##  2 Transfer Function Analysis of the Linear Feedback System

Let $G$ be a strictly proper, asymptotically stable, discrete-time SISO
transfer function with a zero at 1 and no other zeros on the unit circle. Let
$G=N/D$, where the polynomials $N$ and $D$ are coprime, $D$ is monic,
$n\stackrel{{\scriptstyle\triangle}}{{=}}\deg D$, and
$m\stackrel{{\scriptstyle\triangle}}{{=}}\deg N<n.$ Note that $N(1)=0$ and
$D(1)\neq 0,$ and thus $G(1)=0.$ Furthermore, for all
$\theta\in(-\pi,\pi]\backslash\\{0\\},$ $G(e^{\jmath\theta})\neq 0.$

For all $\alpha\in{\mathbb{R}},$ the closed-loop transfer function from $v$ to
$y$ of the linear feedback system in Fig. 1 is given by

$G_{yv,\alpha}({\bf q})\stackrel{{\scriptstyle\triangle}}{{=}}\frac{\alpha
G({\bf q})}{1-\alpha G({\bf q})}=\frac{\alpha N({\bf q})}{p_{\alpha}({\bf
q})},$ (1)

where $p_{\alpha}({\bf q})\stackrel{{\scriptstyle\triangle}}{{=}}D({\bf
q})-\alpha N({\bf q})$. The forward shift operator ${\bf q}$ accounts for both
the free and forced response of the linear feedback system in Fig. 1; for
pole-zero analysis, ${\bf q}$ is replaced by the Z-transform complex variable
${\bf z}$. Note that $p_{0}=D,$ and thus $\operatorname{spr}(p_{0})<1.$
However, for all $|\alpha|$ sufficiently large, it follows from the root locus
asymptote rule that $p_{\alpha}$ has at least $n-m$ roots outside the closed
unit disk, and thus $\operatorname{spr}(p_{\alpha})>1.$ Note that, since Fig.
1 has no sign change in the loop, the root locus parameter $\alpha$ plays the
role of $-k$ in the standard root locus. The following result is immediate.

$+$$\alpha G({\bf q})$$v$$\nu$$y$

Figure 1: Discrete-time linear feedback system with input $v$ and output $y$.

###### Proposition 2.1

Let $\alpha\in{\mathbb{R}}\backslash\\{0\\}$ and
$\theta\in(-\pi,\pi]\backslash\\{0\\}$. Then $p_{\alpha}(e^{\jmath\theta})=0$
if and only if $\alpha=1/G(e^{\jmath\theta}).$

Proposition 2.1 implies that, if $\theta\in(-\pi,\pi]\backslash\\{0\\}$ and
$G(e^{\jmath\theta})$ is real, then $e^{\jmath\theta}$ is a pole of
$G_{yv,1/G(e^{\jmath\theta})}$ and thus an element of either the 0-deg or
180-deg root locus of $G_{yv,\alpha}$. Now, define

$\displaystyle\Theta_{\rm n}$
$\displaystyle\stackrel{{\scriptstyle\triangle}}{{=}}\\{\theta\in(-\pi,\pi]\backslash\\{0\\}\colon
G(e^{\jmath\theta})\in(-\infty,0)\\},$ (2) $\displaystyle\Theta_{\rm p}$
$\displaystyle\stackrel{{\scriptstyle\triangle}}{{=}}\\{\theta\in(-\pi,\pi]\backslash\\{0\\}\colon
G(e^{\jmath\theta})\in(0,\infty)\\},$ (3)

so that $\Theta_{\rm n}\ \cup\ \Theta_{\rm
p}=\\{\theta\in(-\pi,\pi]\backslash\\{0\\}\colon G(e^{\jmath\theta})\mbox{ is
real}\\}.$ Note that $\Theta_{\rm n}$ is the set of angles at which the
180-deg root locus of $G_{yv,\alpha}$ crosses the unit circle, and
$\Theta_{\rm p}$ is the set of angles at which the 0-deg root locus of
$G_{yv,\alpha}$ crosses the unit circle, which occurs in both cases for
$\alpha=1/G(e^{\jmath\theta}).$ Since the 180-deg and 0-deg root locus plots
of $G_{yv,\alpha}$ have $n-m$ asymptotes as $\alpha\to-\infty$ and
$\alpha\to\infty$, respectively, it follows that

$\displaystyle\operatorname{card}(\Theta_{\rm n})$ $\displaystyle\geq n-m,$
(4) $\displaystyle\operatorname{card}(\Theta_{\rm p})$ $\displaystyle\geq
n-m.$ (5)

Furthermore, in the case where $n-m=1,$ the positive real axis is an asymptote
of the 0-deg root locus plot. Since $G$ has a zero at 1, it follows that two
poles must break in on the positive real axis, which implies that

$\operatorname{card}(\Theta_{\rm p})\geq\min\\{2,n-m\\},$ (6)

as illustrated by the following example.

###### Example 2.2

Let $\alpha\geq 0$ and $G({\bf q})=\tfrac{({\bf q}-0.2\pm\jmath 0.8)({\bf
q}-1)}{{\bf q}^{2}({\bf q}+0.25\pm\jmath 0.7)}=\tfrac{({\bf q}^{2}-0.4{\bf
q}+0.68)({\bf q}-1)}{{\bf q}^{2}({\bf q}^{2}+0.5{\bf q}+0.5525)}.$ Since $n=4$
and $m=3,$ it follows that the 0-deg root locus of the closed-loop system has
one asymptote, as shown in Fig. 2. However, the root locus plot crosses the
unit circle at two points due the pole break-in on the positive real axis, and
thus ${\rm card}(\Theta_{\rm p})=2.$

Figure 2: Example 2.2: 0-deg root locus of the linear feedback system in Fig.
1.

The following result implies that the root locus plot of $G$ intersects the
unit circle at a finite number of points.

###### Proposition 2.3

$\Theta_{\rm p}$ and $\Theta_{\rm n}$ are finite.

Proof. Let $r$ be a positive integer such that $h({\bf
z})\stackrel{{\scriptstyle\triangle}}{{=}}{\bf z}^{r}[D({\bf z})N(1/{\bf
z})-N({\bf z})D(1/{\bf z})]$ is a polynomial. Now, let ${\bf
z}=e^{\jmath\theta}$, where $\theta\in\Theta_{\rm n}\cup\Theta_{\rm p}.$ Since
$G({\bf z})$ is real and $|{\bf z}|=1$, it follows that

$\frac{D({\bf z})}{N({\bf z})}=\overline{\left(\frac{D({\bf z})}{N({\bf
z})}\right)}={\frac{\overline{D}({\bf z})}{\overline{N}({\bf
z})}}=\frac{D(\bar{{\bf z}})}{N(\bar{{\bf z}})}=\frac{D(1/{\bf z})}{N(1/{\bf
z})}.$

Hence, $h({\bf z})=0.$ Since $h$ has a finite number of roots, it follows that
$\Theta_{\rm n}$ and $\Theta_{\rm p}$ are finite. $\square$

Proposition 2.3 implies that we can define

$\displaystyle\alpha_{\rm n}$
$\displaystyle\stackrel{{\scriptstyle\triangle}}{{=}}\max_{\theta\in\Theta_{{\rm
n}}}1/G(e^{\jmath\theta})<0,$ (7) $\displaystyle\alpha_{\rm p}$
$\displaystyle\stackrel{{\scriptstyle\triangle}}{{=}}\min_{\theta\in\Theta_{{\rm
p}}}1/G(e^{\jmath\theta})>0.$ (8)

###### Proposition 2.4

If $\alpha\in(\alpha_{\rm n},\alpha_{\rm p}),$ then $G_{yv,\alpha}$ is
asymptotically stable. Furthermore,

$\displaystyle\operatorname{spr}(p_{\alpha_{\rm
n}})=\operatorname{spr}(p_{\alpha_{\rm p}})=1.$ (9)

Proof. Suppose there exists $\alpha\in(\alpha_{\rm n},0)$ such that
$\operatorname{spr}(p_{\alpha})\geq 1.$ Since $\operatorname{spr}$ is
continuous and $\operatorname{spr}(p_{0})=\operatorname{spr}(D)<1,$ the
intermediate value theorem implies that there exists
$\alpha_{1}\in(\alpha_{\rm n},0)$ such that
$\operatorname{spr}(p_{\alpha_{1}})=1.$ Hence, there exists
$\theta_{1}\in(-\pi,\pi]/\\{0\\}$ such that
$p_{\alpha_{1}}(e^{\jmath\theta_{1}})=0.$ Hence,
$G(e^{\jmath\theta_{1}})=1/\alpha_{1}<0,$ and thus $\theta_{1}\in\Theta_{\rm
n}.$ Therefore, (7) implies that

$\max_{\theta\in\Theta_{{\rm n}}}1/G(e^{\jmath\theta})=\alpha_{\rm
n}<\alpha_{1}=1/G(e^{\jmath\theta_{1}}),$

which is a contradiction. Hence, for all $\alpha\in(\alpha_{\rm n},0],$
$G_{yv,\alpha}$ is asymptotically stable. Similarly, for all
$\alpha\in[0,\alpha_{\rm p}),$ $G_{yv,\alpha}$ is asymptotically stable.
Hence, for all $\alpha\in(\alpha_{\rm n},\alpha_{\rm p}),$ $G_{yv,\alpha}$ is
asymptotically stable.

Next, let $\theta_{\rm n}\in\Theta_{\rm n}$ satisfy $\alpha_{\rm
n}=1/G(e^{\jmath\theta_{\rm n}}).$ Then, Proposition 2.1 implies that
$p_{\alpha_{\rm n}}(e^{\jmath\theta_{\rm n}})=0,$ that is,
$e^{\jmath\theta_{\rm n}}$ is a root of $p_{\alpha_{\rm n}},$ and thus
$\operatorname{spr}(p_{\alpha_{\rm n}})\geq 1.$ Now, suppose that
$\operatorname{spr}(p_{\alpha_{\rm n}})>1.$ Since
$\operatorname{spr}(p_{0})<1$ and $\operatorname{spr}$ is continuous, it
follows that there exists $\alpha_{1}\in(\alpha_{\rm n},0)$ such that
$\operatorname{spr}(p_{\alpha_{1}})=1.$ Since $\alpha_{1}\in(\alpha_{\rm
n},\alpha_{\rm p}),$ it follows that $G_{yv,\alpha_{1}}$ is asymptotically
stable, and thus $\operatorname{spr}(p_{\alpha_{1}})<1,$ which is a
contradiction. Hence, $\operatorname{spr}(p_{\alpha_{\rm n}})=1.$ Similarly,
$\operatorname{spr}(p_{\alpha_{\rm p}})=1.$ $\square$

The following result is an immediate consequence of the root locus asymptote
rule.

###### Lemma 2.5

There exist $\beta_{\rm n}<0$ and $\beta_{\rm p}>0$ such that, for all
$\alpha<\beta_{\rm n}$ and all $\alpha>\beta_{\rm p},$ $p_{\alpha}$ has at
least $n-m$ roots with absolute value greater than 1, all of which are simple.

The following example shows that $\alpha_{\rm p}$ defined by (8) is not
necessarily the supremum of all values of $\alpha$ such that $G_{yv,\alpha}$
is asymptotically stable. In other words, there may exist $\alpha>\alpha_{\rm
p}$ such that $G_{yv,\alpha}$ is asymptotically stable.

###### Example 2.6

Let $G({\bf q})=\frac{({\bf q}-0.05\pm\jmath 0.88)({\bf q}-1)}{{\bf
q}^{2}({\bf q}+0.05\pm\jmath 0.88)}=\frac{({\bf q}^{2}-0.1{\bf q}+0.7769)({\bf
q}-1)}{{\bf q}^{2}({\bf q}^{2}+0.1{\bf q}+0.7769)}.$ Fig. 3a) shows the root
locus for $\alpha>0$ and $\alpha<0$, and Fig. 3b) shows
$\operatorname{spr}(p_{\alpha})$ versus $\alpha$, which indicates that there
exists $\alpha>\alpha_{\rm p}$ such that $\operatorname{spr}(p_{\alpha})<1.$

Figure 3: Example 2.6: Root locus and $\operatorname{spr}(p_{\alpha})$ for the
linear feedback system in Fig. 1. a) shows the root locus, where red
corresponds to $\alpha>0$ (the 0-deg root locus), and blue corresponds to
$\alpha<0$ (the 180-deg root locus). b) shows $\operatorname{spr}(p_{\alpha})$
versus $\alpha$, where $\alpha_{\rm n}$ and $\alpha_{\rm p}$ are indicated.
Note that $\alpha_{\rm p}\approx 0.6$, whereas the closed-loop system is
asymptotically stable for all $\alpha\in(1.05,1.2)$.

##  3 State Space Analysis of the Linear Feedback System

Let $(A,B,C)$ be a minimal realization of $G$ with state
$x_{k}\in{\mathbb{R}}^{n}$ at step $k$ so that $A$ is asymptotically stable.
The linear feedback system in Fig. 1 has the closed-loop dynamics

$\displaystyle x_{k+1}$ $\displaystyle=(A+\alpha BC)x_{k}+\alpha Bv_{k},$ (10)
$\displaystyle y_{k}$ $\displaystyle=Cx_{k}.$ (11)

Note that, for all $\alpha\neq 0,$ $(A+\alpha BC,\alpha B,C)$ is a minimal
realization of $G_{yv,\alpha},$ and thus $p_{\alpha}=D-\alpha N=\chi_{A+\alpha
BC}.$

###### Lemma 3.1

Let $\xi\in{\mathbb{C}}^{n}$ be an eigenvector of $A+\alpha BC.$ Then,
$C\xi\neq 0.$

Proof. Suppose that $C\xi=0.$ Since $\xi\in{\mathbb{C}}^{n}$ is an eigenvector
of $A+\alpha BC,$ it follows that $\xi$ is an eigenvector of $A.$ Hence, for
all $i\in\\{0,1,\ldots,n-1\\}$, $CA^{i}\xi=0.$ Since $(A,C)$ is observable, it
follows that $\xi=0,$ which is a contradiction. $\square$

The following lemma concerns the linear feedback system in Fig. 1 with
$v\equiv 0.$ In this case, (10) and (11) can be written as

$\displaystyle x_{k}$ $\displaystyle=(A+\alpha BC)^{k}x_{0},$ (12)
$\displaystyle y_{k}$ $\displaystyle=Cx_{k}.$ (13)

Recall from Lemma 2.5 that, for all $|\alpha|$ sufficiently large, $A+\alpha
BC$ has at least $n-m$ eigenvalues with absolute value greater than 1, all of
which are simple.

###### Lemma 3.2

Consider the linear feedback system in Fig. 1 with $v\equiv 0.$ Let
$\alpha\in{\mathbb{R}},$ assume there exists a simple eigenvalue
$\lambda\in\operatorname{spec}(A+\alpha BC)$ such that $|\lambda|>1$, let
$\xi\in{\mathbb{C}}^{n}$ be an associated eigenvector, let
${\mathcal{X}}\subset{\mathbb{C}}^{n}$ be the $n-1$-dimensional subspace
spanned by the eigenvectors and generalized eigenvectors associated with the
remaining eigenvalues of $A+\alpha BC,$ and assume that
$x_{0}\notin{\mathcal{X}}.$ Then the following statements hold:

  1. $i)$

For all $k\geq 0$, $x_{k}\notin{\mathcal{X}}.$

  2. $ii)$

$\limsup_{k\to\infty}|y_{k}|=\infty.$

Proof. Since ($A,B$) is controllable, Fact 7.15.10 of [18, p. 599] implies
that $A$ is cyclic. Since, in addition, ($A+\alpha BC,B$) is controllable, it
follows that $A+\alpha BC$ is cyclic. Therefore, each eigenvalue of $A+\alpha
BC$ has exactly one associated eigenvector. Let
$\lambda_{1},\ldots,\lambda_{r}$ be the distinct eigenvalues of $A+\alpha BC,$
for all $j\in\\{1,\ldots,r\\},$ let $n_{j}$ be the algebraic multiplicity of
$\lambda_{j},$ and, for all $i\in\\{1,\ldots,n_{j}\\},$ let $\xi_{j,i}$ be a
generalized eigenvector of $A+\alpha BC$ corresponding to $\lambda_{j}$ such
that $(A+\alpha BC-\lambda_{j}I)^{i-1}\xi_{j,i}\neq 0,$ $(A+\alpha
BC-\lambda_{j}I)^{i}\xi_{j,i}=0,$ and $(\xi_{j,1},\ldots,\xi_{j,n_{j}})$ is a
Jordan chain of $A+\alpha BC$ associated with $\lambda_{j}.$ Note that, for
all $j\in\\{1,\ldots,r\\},$ $\xi_{j,1}$ is an eigenvector associated with
$\lambda_{j}$. Without loss of generality, define
$\lambda_{1}\stackrel{{\scriptstyle\triangle}}{{=}}\lambda$ and
$\xi_{1,1}\stackrel{{\scriptstyle\triangle}}{{=}}\xi.$ Note that, since
$\lambda$ is simple, it follows that $n_{1}=1.$

Next, it follows from the Jordan decomposition and equation (7.8.5) from [19,
p. 594] that

$(A+\alpha BC)^{k}=S\begin{bmatrix}J_{1}^{k}&&\\\ &\ddots&\\\
&&J_{r}^{k}\end{bmatrix}S^{-1},$ (14)

where $S\stackrel{{\scriptstyle\triangle}}{{=}}[\xi_{1}\ \cdots\
\xi_{r}]\in{\mathbb{C}}^{n\times n},$ for all $j\in\\{1,\ldots,r\\},$
$\xi_{j}\stackrel{{\scriptstyle\triangle}}{{=}}[\xi_{j,1}\ \cdots\
\xi_{j,n_{j}}]\in{\mathbb{C}}^{n\times n_{j}},$
$J_{j}\in{\mathbb{C}}^{n_{j}\times n_{j}}$ is the Jordan block associated with
the eigenvalue $\lambda_{j}$ of $A+\alpha BC,$ and

$J_{j}^{k}=\begin{bmatrix}\lambda_{j}^{k}&\binom{k}{1}\lambda_{j}^{k-1}&\binom{k}{2}\lambda_{j}^{k-2}&\cdots&\cdots&\binom{k}{n_{j}-1}\lambda_{j}^{k-n_{j}+1}\\\
&\lambda_{j}^{k}&\binom{k}{1}\lambda_{j}^{k-1}&\cdots&\cdots&\binom{k}{n_{j}-2}\lambda_{j}^{k-n_{j}+2}\\\
&&\ddots&\ddots&\vdots&\vdots\\\ &&&\ddots&\ddots&\vdots\\\
&&&&\lambda_{j}^{k}&\binom{k}{1}\lambda_{j}^{k-1}\\\
&&&&&\lambda_{j}^{k}\end{bmatrix},$

(15)

where, for all $k<i,$ $\binom{k}{i}\stackrel{{\scriptstyle\triangle}}{{=}}0.$
Since $S$ is invertible, the $n$ linearly independent generalized eigenvectors
$\xi,\xi_{2,1},\ldots,\xi_{2,n_{2}},\ldots,\xi_{r,1},\ldots,\xi_{r,n_{r}}$
comprise a basis of ${\mathbb{C}}^{n}$. Therefore, it follows that, for all
$j\in\\{1,\ldots,r\\}$ and $i\in\\{1,\ldots,n_{j}\\},$ there exists
$\beta_{j,i}\in{\mathbb{C}}$ such that

$\displaystyle x_{0}=\sum_{j=1}^{r}\sum_{i=1}^{n_{j}}\beta_{j,i}\xi_{j,i}$
$\displaystyle=\beta_{1,1}\xi+\sum_{j=2}^{r}\sum_{i=1}^{n_{j}}\beta_{j,i}\xi_{j,i}$
$\displaystyle=S\begin{bmatrix}\beta_{1}^{\rm T}&\dots&\beta_{r}^{\rm
T}\end{bmatrix}^{\rm T},$ (16)

where $\beta_{j}\stackrel{{\scriptstyle\triangle}}{{=}}[\beta_{j,1}\ \cdots\
\beta_{j,n_{j}}]^{\rm T}\in{\mathbb{C}}^{n_{j}\times 1}.$ It thus follows from
(12), (14), (15), and (16) that, for all $k\geq 0,$

$\displaystyle x_{k}$ $\displaystyle=S\begin{bmatrix}J_{1}^{k}&&\\\
&\ddots&\\\ &&J_{r}^{k}\end{bmatrix}\begin{bmatrix}\beta_{1}\\\ \vdots\\\
\beta_{r}\end{bmatrix}=S\begin{bmatrix}J_{1}^{k}\beta_{1}\\\ \vdots\\\
J_{r}^{k}\beta_{r}\end{bmatrix}$
$\displaystyle=\sum_{j=1}^{r}\sum_{i=1}^{n_{j}}\gamma_{j,i}(k)\lambda_{j}^{k}\xi_{j,i},$
(17)

where, for all $j\in\\{1,\ldots,r\\}$, $i\in\\{1,\ldots,n_{j}\\},$ and $k\geq
0,$

$\gamma_{j,i}(k)\stackrel{{\scriptstyle\triangle}}{{=}}\beta_{j,i}+\sum_{l=2}^{n_{j}-i+1}\binom{k}{l-1}\lambda_{j}^{1-l}\beta_{j,l+i-1}.$
(18)

In particular,

$\gamma_{1,1}(k)=\beta_{1,1}.$ (19)

Hence,

$\displaystyle x_{k}$
$\displaystyle=\beta_{1,1}\lambda^{k}\xi+\sum_{j=2}^{r}\sum_{i=1}^{n_{j}}\gamma_{j,i}(k)\lambda_{j}^{k}\xi_{j,i},$
(20)

and thus it follows from (13) that

$y_{k}=p_{1}(k)\lambda^{k}+\sum_{j=2}^{r}p_{j}(k)\lambda_{j}^{k},$ (21)

where

$\displaystyle p_{1}(k)$
$\displaystyle\stackrel{{\scriptstyle\triangle}}{{=}}\beta_{1,1}C\xi,$ (22)

and, for all $j=2,\ldots,r,$

$p_{j}(k)\stackrel{{\scriptstyle\triangle}}{{=}}\sum_{i=1}^{n_{j}}\left(\beta_{j,i}+\sum_{l=2}^{n_{j}-i+1}\binom{k}{l-1}\beta_{j,l+i-1}\lambda_{j}^{1-l}\right)C\xi_{j,i}.$

(23)

Note that $p_{1}(k),\ldots,p_{r}(k)$ are polynomials in $k$ with complex
coefficients. Now, the $n-1$-dimensional subspace spanned by the eigenvectors
and generalized eigenvectors associated with the eigenvalues
$\lambda_{2},\ldots,\lambda_{r}$ of $A+\alpha BC$ is given by

${\mathcal{X}}={\mathcal{R}}([\xi_{2,1}\ \cdots\ \xi_{2,n_{2}}\ \cdots\
\xi_{r,1}\cdots\xi_{r,n_{r}}]).$ (24)

It follows from (16) that, since $x_{0}\notin{\mathcal{X}},$ $\beta_{1,1}\neq
0,$ and thus it follows from (20) that, for all $k\geq 0,$
$x_{k}\notin{\mathcal{X}},$ which proves i). Next, Lemma 3.1 implies that
$C\xi\neq 0.$ Hence, $p_{1}$ is not the zero polynomial. Finally, since
$p_{1}$ is not the zero polynomial and $|\lambda|>1,$ Proposition A.3 implies
that $\limsup_{k\to\infty}\left|y_{k}\right|=\infty,$ which proves ii).
$\square$

The following alternative characterization of ${\mathcal{X}}$ is worth noting.

###### Proposition 3.3

Let $\alpha\in{\mathbb{R}},$ assume there exists a simple eigenvalue
$\lambda\in\operatorname{spec}(A+\alpha BC)$ such that $|\lambda|>1$, define
${\mathcal{X}}$ as in Lemma 3.2, and define the polynomial
$q(z)\stackrel{{\scriptstyle\triangle}}{{=}}\chi_{A+\alpha
BC}(z)/(z-\lambda).$ Then ${\mathcal{X}}={\mathcal{N}}(q(A+\alpha BC)).$

The following is a corollary of Lemma 3.2.

###### Corollary 3.4

Consider the linear feedback system in Fig. 1 with $v\equiv 0.$ Let
$\alpha\in{\mathbb{R}},$ assume there exists a simple eigenvalue
$\lambda\in\operatorname{spec}(A+\alpha BC)$ such that $|\lambda|>1$, let
$\xi\in{\mathbb{C}}^{n}$ be an associated eigenvector, let
${\mathcal{X}}\subset{\mathbb{C}}^{n}$ be the $n-1$-dimensional subspace
spanned by the eigenvectors and generalized eigenvectors associated with the
remaining eigenvalues of $A+\alpha BC$, let $k_{0}\geq 0,$ assume that
$x_{k_{0}}\notin{\mathcal{X}},$ and let $M>|y_{k_{0}}|$. Then, the following
statements hold:

  1. $i)$

For all $k\geq k_{0}$, $x_{k}\notin{\mathcal{X}}.$

  2. $ii)$

There exists $k_{1}>k_{0}$ such that the following statements hold:
$a$) For all $k\in\\{k_{0},\ldots,k_{1}-1\\},$ $|y_{k}|<M$.
$b$) $|y_{k_{1}}|\geq M.$

##  4 Analysis of the Lur’e System

We now consider the discrete-time Lur’e system in Fig. 4, which has the
closed-loop dynamics

$\displaystyle x_{k+1}$ $\displaystyle=Ax_{k}+\alpha B{\rm sat}_{1}(y_{k}),$
(25) $\displaystyle y_{k}$ $\displaystyle=Cx_{k}.$ (26)

$+$$\alpha G({\bf q})$${\rm sat}_{1}$$v\equiv 0$$\nu$$y$

Figure 4: Discrete-time Lur’e system.

###### Proposition 4.1

Let $y$ be the output of the discrete-time Lur’e system in Fig. 4. Then, $y$
is bounded.

Proof: Since $G$ is asymptotically stable, it is bounded-input, bounded-output
stable. Since the input of $G$ is bounded, it follows that $y$ is bounded.
$\square$

Let $S_{k}$ denote the dynamics of the discrete-time Lur’e system at step
$k\geq 0,$ and consider the transition diagram in Fig. 5, where, for all
$k\geq 0$ such that $y_{k}\geq 1,$ $S_{k}={\mathcal{S}}_{1}$, and thus
$x_{k+1}=Ax_{k}+\alpha B;$ for all $k\geq 0$ such that $|y_{k}|<1,$
$S_{k}={\mathcal{S}}_{2}$, and thus $x_{k+1}=(A+\alpha BC)x_{k};$ and, for all
$k\geq 0$ such that $y_{k}\leq-1,$ $S_{k}={\mathcal{S}}_{3}$, and thus
$x_{k+1}=Ax_{k}-\alpha B.$

${\mathcal{S}}_{2}$$+$$\alpha G({\bf
q})$$v=0$$\nu$$y$${\mathcal{S}}_{1}$$\alpha G({\bf
q})$$\nu=1$$y$${\mathcal{S}}_{3}$$\alpha G({\bf q})$$\nu=-1$$y$$y_{k}\geq
1$$y_{k}\in(-1,1)$$y_{k}\leq-1$$y_{k}\in(-1,1)$$y_{k}\geq
1$$y_{k}\leq-1$$y_{k}\geq 1$$y_{k}\leq-1$

Figure 5: Transition diagram for the discrete-time Lur’e system.

Henceforth, the terminology $\lim_{k\to\infty}y_{k}$ does not exist means that
there does not exist a real number $\overline{y}$ such that
$\lim_{k\to\infty}y_{k}=\overline{y}.$ For example,
$\lim_{k\to\infty}(-1)^{k}$, $\lim_{k\to\infty}k,$ and
$\lim_{k\to\infty}(-1)^{k}k$ do not exist. The following lemma is needed.

###### Lemma 4.2

Let $\alpha\in{\mathbb{R}}.$ Then the following statements hold:

  1. $i)$

$1\notin\operatorname{spec}(A+\alpha BC)$.

  2. $ii)$

If $\lim_{k\to\infty}(A+\alpha BC)^{k}$ exists, then the limit is zero.

  3. $iii)$

If, for $y_{k}$ given by (10) and (11) with $v_{k}\equiv 0,$
$\lim_{k\to\infty}y_{k}$ exists, then the limit is zero.

Proof: To prove i), first consider the case where $\alpha=0.$ Since $A$ is
asymptotically stable, it follows that $1\notin\operatorname{spec}(A)$. Now,
let $\alpha\neq 0.$ Then $(A+\alpha BC,\alpha B,C)$ is a minimal realization
of $G_{yv,\alpha},$ and thus $\chi_{A+\alpha BC}=p_{\alpha}=D-\alpha N.$ Since
$N$ and $D$ are coprime and $G(1)=0$, it follows that $N(1)=0$ and $D(1)\neq
0.$ Hence, $\chi_{A+\alpha BC}(1)=D(1)-\alpha N(1)=D(1)\neq 0.$ Therefore,
$1\notin\operatorname{spec}(A+\alpha BC)$.

ii) and iii) follow from Lemma A.4. $\square$

###### Lemma 4.3

Let $y$ be the output of the discrete-time Lur’e system in Fig. 4 with
$\alpha\in{\mathbb{R}}\backslash\\{0\\}$. Consider the following statements:
$a$) $\lim_{k\to\infty}y_{k}$ exists, and $\lim_{k\to\infty}y_{k}=0.$
$b$) $\lim_{k\to\infty}y_{k}$ exists, and $\lim_{k\to\infty}y_{k}\in(-1,1).$
$c$) $\lim_{k\to\infty}y_{k}$ exists, and $\lim_{k\to\infty}y_{k}\in[-1,1].$
$d$) $\lim_{k\to\infty}y_{k}$ exists.
Then $a$) $\Longleftrightarrow$ $b$) $\Longleftrightarrow$ $c$)
$\Longleftrightarrow$ $d$).

Proof: Note that $a$) $\Longrightarrow$ $b$) $\Longrightarrow$ $c$)
$\Longrightarrow$ $d$).

To prove $b)$ $\Longrightarrow$ $a),$ note that
$\lim_{k\to\infty}y_{k}\in(-1,1)$ implies that there exists $k_{1}\geq n$ such
that, for all $k\geq k_{1}$, $|y_{k}|<1$ and $S_{k}={\mathcal{S}}_{2}.$ Hence,
for all $k\geq k_{1}$ $y_{k}$ is given by (12) and (13) with $x_{0}$ replaced
by $x_{k_{1}}$. Since $\lim_{k\to\infty}y_{k}=\lim_{k\to\infty}C(A+\alpha
BC)^{k-k_{1}}x_{k_{1}}$ exists, Lemma 4.2 implies that
$\lim_{k\to\infty}y_{k}=0.$

To prove $c)$ $\Longrightarrow$ $b),$ consider the case where
$\lim_{k\to\infty}y_{k}=1,$ so that $\lim_{k\to\infty}\nu_{k}=1.$ Next, note
that it follows from input-to-state stability for linear time-invariant
discrete-time systems (see [20, Example 3.4]) that

$\lim_{k\to\infty}\sum_{i=0}^{k-1}CA^{k-1-i}B(\nu_{i}-1)=0.$ (27)

Since $\lim_{k\to\infty}CA^{k}x_{0}=0$ and $C(I-A)^{-1}B=G(1)=0,$ it follows
from (27) that

$\displaystyle 1$ $\displaystyle=\lim_{k\to\infty}y_{k}$
$\displaystyle=\lim_{k\to\infty}CA^{k}x_{0}+\alpha\lim_{k\to\infty}\sum_{i=0}^{k-1}CA^{k-1-i}B\nu_{i}$
$\displaystyle=\alpha\lim_{k\to\infty}\sum_{i=0}^{k-1}CA^{k-1-i}B$
$\displaystyle\quad+\alpha\lim_{k\to\infty}\sum_{i=0}^{k-1}CA^{k-1-i}B(\nu_{i}-1)$
$\displaystyle=\alpha\sum_{i=0}^{\infty}CA^{i}B=\alpha\ C(I-A)^{-1}B$
$\displaystyle=0,$

which is a contradiction. Hence, $\lim_{k\to\infty}y_{k}<1.$ A similar
argument implies that $\lim_{k\to\infty}y_{k}>-1.$ Hence,
$\lim_{k\to\infty}y_{k}\in(-1,1).$

To prove $d$) $\Longrightarrow$ $c$), suppose that $\lim_{k\to\infty}y_{k}>1.$
Then there exists $k_{1}\geq n$ such that, for all $k\geq k_{1}$, $y_{k}\geq
1,$ and thus $S_{k}={\mathcal{S}}_{1}.$ Hence, for all $k\geq k_{1},$ $y_{k}$
is given by

$\displaystyle y_{k}$
$\displaystyle=CA^{k-k_{1}}x_{k_{1}}+\alpha\sum_{i=0}^{k-k_{1}-1}CA^{k-k_{1}-1-i}B$
$\displaystyle=CA^{k-k_{1}}x_{k_{1}}+\alpha\sum_{i=0}^{k-k_{1}-1}CA^{i}B.$

Since $\lim_{k\to\infty}A^{k}=0$ and $C(I-A)^{-1}B=G(1)=0,$ it follows that

$\displaystyle\lim_{k\to\infty}y_{k}$
$\displaystyle=\lim_{k\to\infty}CA^{k-k_{1}}x_{k_{1}}+\lim_{k\to\infty}\alpha\sum_{i=0}^{k-k_{1}-1}CA^{i}B$
$\displaystyle=\alpha\sum_{i=0}^{\infty}CA^{i}B=\alpha C(I-A)^{-1}B=0.$

Hence, there exists $k_{2}>k_{1}$ such that $y_{k_{2}}<1,$ which is a
contradiction. Therefore, $\lim_{k\to\infty}y_{k}\leq 1.$ Similarly,
$\lim_{k\to\infty}y_{k}\geq-1.$ Hence, $\lim_{k\to\infty}y_{k}\in[-1,1].$
$\square$

###### Lemma 4.4

Consider the discrete-time Lur’e system in Fig. 5. Let
$\alpha\in{\mathbb{R}},$ assume there exists a simple eigenvalue
$\lambda\in\operatorname{spec}(A+\alpha BC)$ such that $|\lambda|>1$, let
$\xi\in{\mathbb{C}}^{n}$ be an associated eigenvector, let
${\mathcal{X}}\subset{\mathbb{C}}^{n}$ be the $n-1$-dimensional subspace
spanned by the eigenvectors and generalized eigenvectors associated with the
remaining eigenvalues of $A+\alpha BC,$ let $k_{0}\geq 0$, and assume that
$x_{k_{0}}\notin{\mathcal{X}}$ and $S_{k_{0}}={\mathcal{S}}_{2}$. Then there
exists $k_{1}>k_{0}$ such that, for all
$k\in\\{k_{0},k_{0}+1,\ldots,k_{1}-1\\},$ $S_{k}={\mathcal{S}}_{2}$ and
$x_{k}\notin{\mathcal{X}},$ and such that $x_{k_{1}}\notin{\mathcal{X}}$ and
$S_{k_{1}}\neq{\mathcal{S}}_{2}.$

Proof: For all $k\geq k_{0}$ such that $S_{k}={\mathcal{S}}_{2},$ the system
dynamics are given by the discrete-time linear feedback system in Fig. 1 with
$v\equiv 0.$ Hence, the result follows from Corollary 3.4 with $M=1.$
$\square$

The following result gives sufficient conditions under which the response of
the discrete-time Lur’e system is bounded and nonconvergent.

###### Theorem 4.5

Consider the discrete-time Lur’e system in Fig. 5. Let
$\alpha\in{\mathbb{R}},$ assume there exists a simple eigenvalue
$\lambda\in\operatorname{spec}(A+\alpha BC)$ such that $|\lambda|>1$, let
$\xi\in{\mathbb{C}}^{n}$ be an associated eigenvector, and let
${\mathcal{X}}\subset{\mathbb{C}}^{n}$ be the $n-1$-dimensional subspace
spanned by the eigenvectors and generalized eigenvectors associated with the
remaining eigenvalues of $A+\alpha BC$. Furthermore, assume that, if
$S_{0}={\mathcal{S}}_{2},$ then $x_{0}\notin{\mathcal{X}}$, and, for all
$k\geq 0$ such that $S_{k}\neq{\mathcal{S}}_{2}$ and
$S_{k+1}={\mathcal{S}}_{2}$, it follows that $x_{k+1}\notin{\mathcal{X}}.$
Then $y$ is bounded and $\lim_{k\to\infty}y_{k}$ does not exist.

Proof: Proposition 4.1 implies that $y$ is bounded. Suppose that
$\lim_{k\to\infty}y_{k}$ exists, and thus Lemma 4.3 implies that
$\lim_{k\to\infty}y_{k}=0.$ Now, let $k_{0}\geq 0$ be the smallest nonnegative
integer such that, for all $k\geq k_{0},$ $|y_{k}|<1$ and thus
$S_{k}={\mathcal{S}}_{2}.$ In the case where $k_{0}=0,$ it follows that
$x_{0}\notin{\mathcal{X}},$ and thus Lemma 4.4 implies that there exists
$k>k_{0}$ such that $|y_{k}|>1.$ Alternatively, in the case where $k_{0}>0$,
note that ${\mathcal{S}}_{k_{0}-1}\neq{\mathcal{S}}_{2}$ and
$x_{k_{0}}\notin{\mathcal{X}},$ and thus Lemma 4.4 implies that there exists
$k>k_{0}$ such that $|y_{k}|>1.$ In both cases, the existence of $k>k_{0}$
such that $|y_{k}|>1$ contradicts the fact that, for all $k\geq k_{0},$
$|y_{k}|<1.$ Therefore, $\lim_{k\to\infty}y_{k}$ does not exist. $\square$

The following result shows that, for almost all initial conditions $x_{0},$
the hypotheses of Theorem 4.5 are satisfied.

###### Theorem 4.6

Consider the discrete-time Lur’e system in Fig. 5. Let
$\alpha\in{\mathbb{R}},$ assume that $A$ is nonsingular, assume there exists a
simple eigenvalue $\lambda\in\operatorname{spec}(A+\alpha BC)$ such that
$|\lambda|>1$, let $\xi\in{\mathbb{C}}^{n}$ be an associated eigenvector, and
let ${\mathcal{X}}\subset{\mathbb{C}}^{n}$ be the $n-1$-dimensional subspace
spanned by the eigenvectors and generalized eigenvectors associated with the
remaining eigenvalues of $A+\alpha BC$. Then, for almost all
$x_{0}\in{\mathbb{R}}^{n},$ $y$ is bounded and $\lim_{k\to\infty}y_{k}$ does
not exist.

Proof: Define $f\colon{\mathbb{R}}^{n}\to{\mathbb{R}}^{n}$ by $f(x)=Ax+\alpha
B{\rm sat}(Cx).$ Letting ${\mathcal{E}}$ denote a proper affine subspace of
${\mathbb{R}}^{n},$ it follows that

$\displaystyle f^{-1}({\mathcal{E}})$ $\displaystyle\subset
A^{-1}({\mathcal{E}}-\alpha B)\cup A^{-1}({\mathcal{E}}+\alpha B)$
$\displaystyle\quad\cup(A+\alpha BC)^{-1}({\mathcal{E}}).$

Hence, the inverse image of every subset of the union of a finite number of
proper affine subspaces of ${\mathbb{R}}^{n}$ is a subset of the union of a
finite number of proper affine subspaces. In particular,
$f^{-1}({\mathcal{X}})$ has measure zero. Now, for all $k\geq 1,$ define
$f^{-k-1}({\mathcal{X}})\stackrel{{\scriptstyle\triangle}}{{=}}f^{-1}(f^{-k}({\mathcal{X}})).$
By induction, it follows that, for all $k\geq 0,$ $f^{-k}({\mathcal{X}})$ a
subset of the union of a finite number of proper affine subspaces and thus has
measure zero. Hence, $\cup_{k\geq 1}f^{-k}({\mathcal{X}})$ is a countable
union of sets with measure zero, and thus has measure zero. Therefore, for all
$x_{0}\notin\cup_{k\geq 1}f^{-k}({\mathcal{X}}),$ it follows that, for all
$k\geq 0,$ $x_{k}\notin{\mathcal{X}}.$ $\square$

##  5 Numerical Example

A key assumption in Theorem 4.5 is the requirement that, if, at step $k,$ the
system changes either from ${\mathcal{S}}_{1}$ to ${\mathcal{S}}_{2}$ or from
${\mathcal{S}}_{3}$ to ${\mathcal{S}}_{2},$ then the state $x_{k+1}$ is not
contained in the $n-1$-dimensional subspace ${\mathcal{X}}.$ Under this
assumption, the state of the closed-loop system includes a component from the
unstable subspace of the closed-loop dynamics. Consequently, the norm of the
state ultimately increases, and thus the system eventually leaves
${\mathcal{S}}_{2}.$ Theorem (4.6) implies that this assumption is satisfied
for almost all initial conditions, as demonstrated by the following example,
where the response converges and thus is not self-excited. However, under a
small perturbation of the initial condition, the response is oscillatory,
which illustrates the generic nature of the assumption on the initial
condition.

###### Example 5.1

Let $G({\bf q})=\tfrac{{\bf q}-1}{{\bf q}^{2}-{\bf q}+0.5},$ so that
$\Theta_{\rm n}=\\{\pi\\}$ rad, $\Theta_{\rm p}=\\{\pm{\rm acos}(3/4)\\}$ rad,
$\alpha_{\rm n}=-1.25,$ and $\alpha_{\rm p}=0.5.$ The root locus of the
closed-loop linear system is shown in Fig. 6, along with $\operatorname{spr}$
for a range of $\alpha.$ Note that, for all $\alpha>\alpha_{\rm p},$ both
eigenvalues of $A+\alpha BC$ are unstable, whereas, for all
$\alpha<\alpha_{\rm n},$ one eigenvalue of $A+\alpha BC$ is unstable and one
eigenvalue of $A+\alpha BC$ is asymptotically stable.

Figure 6: Example 5.1. Root locus and $\operatorname{spr}(p_{\alpha})$ for the
closed-loop linear system. a) shows the root locus, where red corresponds to
$\alpha>0$ (the 0-deg root locus), and blue corresponds to $\alpha<0$ (the
180-deg root locus). b) shows $\operatorname{spr}(p_{\alpha})$ versus
$\alpha$, where $\alpha_{\rm n}$ and $\alpha_{\rm p}$ are indicated.

Let the minimal realization of $G$ be given by

$A=\begin{bmatrix}1&-0.5\\\ 1&0\end{bmatrix},\ B=\begin{bmatrix}1\\\
0\end{bmatrix},\ C=\begin{bmatrix}1&-1\end{bmatrix},$

and let $\alpha=-2.5<\alpha_{\rm n}.$ Then $\operatorname{spec}(A+\alpha
BC)=\\{\lambda,\lambda_{2}\\},$ where
$\lambda=-0.75-0.25\sqrt{41}\approx-2.35$ and
$\lambda_{2}=-0.75+0.25\sqrt{41}\approx 0.85.$ Hence, $|\lambda|>1$ and
$|\lambda_{2}|<1.$ Furthermore, let $\xi=[-0.75-0.25\sqrt{41}\ \ \ 1]^{\rm T}$
and $\xi_{2}=[-0.75+0.25\sqrt{41}\ \ \ 1]^{\rm T}$ be eigenvectors of
$A+\alpha BC$ associated with $\lambda$ and $\lambda_{2}$, respectively.
Hence, ${\mathcal{X}}=\\{\alpha\xi_{2}\colon\alpha\in{\mathbb{R}}\\}.$ Now,
let $\Psi\in{\mathbb{C}}^{2}$ satisfy $||\Psi||=1$ and
${\mathcal{X}}=\\{\Psi\\}^{\perp},$ that is, $\Psi^{*}\xi_{2}=0,$ and define
the projector $P\stackrel{{\scriptstyle\triangle}}{{=}}\Psi\Psi^{*}.$ Note
that, for all $x\in{\mathbb{R}}^{2},$ $x\notin{\mathcal{X}}$ if and only if
$||Px||\neq 0.$

Let $x_{0}$ be given by

$x_{0}=[-5.5-6.5\sqrt{41}\ \ \ -51]^{\rm T},$ (28)

such that $y_{0}=Cx_{0}\approx 3.88>1.$ Then, it follows from the system
dynamics with exact symbolic computation that

$\displaystyle x_{1}$ $\displaystyle=Ax_{0}+\alpha B=[17.5-6.5\sqrt{41}\ \ \ \
-5.5-6.5\sqrt{41}]^{\rm T},$ $\displaystyle x_{2}$
$\displaystyle=Ax_{1}+\alpha B=[17.75-3.25\sqrt{41}\ \ \ \
17.5-6.5\sqrt{41}]^{\rm T},$ $\displaystyle x_{3}$
$\displaystyle=Ax_{2}+\alpha B=[6.5\ \ \ \ 17.75-3.25\sqrt{41}]^{\rm T},$
$\displaystyle x_{4}$ $\displaystyle=Ax_{3}+\alpha B=[-4.875+1.625\sqrt{41}\ \
\ \ 6.5]^{\rm T}.$

Hence, $y_{1}=Cx_{1}=23,$ $y_{2}=Cx_{2}=0.25+3.25\sqrt{41}\approx 21.06>1,$
$y_{3}=Cx_{3}=-11.25+3.25\sqrt{41}\approx 9.56>1,$
$y_{4}=Cx_{4}=-11.375+1.625\sqrt{41}\approx-0.97\in(-1,1),$ and thus
$S_{0}=S_{1}=S_{2}=S_{3}={\mathcal{S}}_{1},$ and $S_{4}={\mathcal{S}}_{2}.$
Note that $x_{4}=6.5\xi_{2}.$ Hence, for all $k\geq 4,$

$\displaystyle C(A+\alpha BC)^{k-4}x_{4}$ $\displaystyle=6.5C(A+\alpha
BC)^{k-4}\xi_{2}$ $\displaystyle=6.5\lambda_{2}^{k-4}C\xi_{2}.$ (29)

Since $|Cx_{4}|=6.5|C\xi_{2}|=|-11.375+1.6250\sqrt{41}|\approx 0.97<1,$ it
follows from (29) that, for all $k\geq 4,$

$\displaystyle|C$ $\displaystyle(A+\alpha
BC)^{k-4}x_{4}|=6.5|\lambda_{2}|^{k-4}|C\xi_{2}|$
$\displaystyle=|-11.375+1.6250\sqrt{41}||-0.75-0.25\sqrt{41}|^{k-4}$
$\displaystyle\approx 0.97\ (0.85)^{k-4}<1,$

and thus $S_{k}={\mathcal{S}}_{2}.$ Hence,

$\displaystyle\lim_{k\to\infty}y_{k}$
$\displaystyle=\lim_{k\to\infty}C(A+\alpha BC)^{k-4}x_{4}$ $\displaystyle=6.5\
C\xi_{2}\lim_{k\to\infty}\lambda_{2}^{k-4}=0,$

and, for all $k\geq 4,$

$\displaystyle||Px_{k}||$ $\displaystyle=||P(A+\alpha BC)^{k-4}x_{4}||$
$\displaystyle=6.5\ ||\lambda_{2}||^{k-4}||P\xi_{2}||=0.$

Therefore, $||Px_{4}||=0,$ which implies that $x_{4}\in{\mathcal{X}},$ and
thus the assumptions of Theorem 4.5 are not satisfied. In this case, the
output converges and thus the system does not have self-excited oscillations,
as shown in Fig. 7.

Next, let $x_{0}$ be given by

$x_{0}=[-5.5-6.5\sqrt{41}\ \ \ -51+\varepsilon]^{\rm T},$ (30)

where $\varepsilon\stackrel{{\scriptstyle\triangle}}{{=}}10^{-12},$ which
represents a small perturbation of (28), such that $y_{0}=Cx_{0}\approx
3.88>1.$ With the initial condition (30), it follows that

$\displaystyle x_{1}$ $\displaystyle=[17.5-6.5\sqrt{41}-0.5\varepsilon\ \ \ \
-5.5-6.5\sqrt{41}]^{\rm T},$ $\displaystyle x_{2}$
$\displaystyle=[17.75-3.25\sqrt{41}-0.5\varepsilon$ $\displaystyle\hskip
64.00012pt17.5-6.5\sqrt{41}-0.5\varepsilon]^{\rm T},$ $\displaystyle x_{3}$
$\displaystyle=[6.5-0.25\varepsilon\ \ \ \
17.75-3.25\sqrt{41}-0.5\varepsilon]^{\rm T},$ $\displaystyle x_{4}$
$\displaystyle=[-4.875+1.625\sqrt{41}\ \ \ \ 6.5-0.25\varepsilon]^{\rm T}.$

Hence, $y_{1}=Cx_{1}=23-0.5\varepsilon,$
$y_{2}=Cx_{2}=0.25+3.25\sqrt{41}\approx 21.06>1,$
$y_{3}=Cx_{3}=-11.25+3.25\sqrt{41}+0.25\varepsilon\approx 9.56>1,$
$y_{4}=Cx_{4}=-11.375+1.625\sqrt{41}+0.25\varepsilon\approx-0.97\in(-1,1),$
and thus $S_{0}=S_{1}=S_{2}=S_{3}={\mathcal{S}}_{1},$ and
$S_{4}={\mathcal{S}}_{2}.$ Defining
$\kappa_{1}\stackrel{{\scriptstyle\triangle}}{{=}}0.5-3\sqrt{41}/82$ and
$\kappa_{2}\stackrel{{\scriptstyle\triangle}}{{=}}0.5+3\sqrt{41}/82,$ it
follows that

$\displaystyle\kappa_{1}\xi+\kappa_{2}\xi_{2}$
$\displaystyle=\begin{bmatrix}-\sqrt{41}/8+9\sqrt{41}/328\\\
0.5-3\sqrt{41}/82\end{bmatrix}$
$\displaystyle+\begin{bmatrix}\sqrt{41}/8-9\sqrt{41}/328\\\
0.5+3\sqrt{41}/82\end{bmatrix}=\begin{bmatrix}0\\\ 1\end{bmatrix}.$

Then, $x_{4}=6.5\xi_{2}+[0\ \ \
-0.25\varepsilon]=(-0.25\kappa_{2}\varepsilon)\xi+(6.5-0.25\kappa_{2}\varepsilon)\xi_{2}=\eta_{1}\xi+\eta_{2}\xi_{2},$
where $\eta_{1},\eta_{2}\neq 0.$ Hence, for all $k\geq 4,$

$\displaystyle C(A+\alpha BC)^{k-4}x_{4}$ $\displaystyle=\eta_{1}C(A+\alpha
BC)^{k-4}\xi$ $\displaystyle+\eta_{2}C(A+\alpha BC)^{k-4}\xi_{2}$
$\displaystyle=\eta_{1}C\xi\lambda^{k-4}+\eta_{2}C\xi_{2}\lambda_{2}^{k-4}.$
(31)

Since $\eta_{1}C\xi\neq 0,$ $\eta_{2}C\xi_{2}\neq 0,$
$|\lambda|>|\lambda_{2}|,$ and $|\lambda|>1,$ it follows from (31) and
Proposition A.3 that

$\displaystyle\limsup_{k\to\infty}$ $\displaystyle\left|C(A+\alpha
BC)^{k-4}x_{4}\right|$
$\displaystyle=\limsup_{k\to\infty}\left|\eta_{1}C\xi\lambda^{k-4}+\eta_{2}C\xi_{2}\lambda_{2}^{k-4}\right|=\infty.$

Hence, there exists $k>4$ such that $S_{k}\neq{\mathcal{S}}_{2}.$ Furthermore,
note that

$||Px_{4}||=||\eta_{1}P\xi+\eta_{2}P\xi_{2}||=||\eta_{1}P\xi||\neq 0,$

which implies that $x_{4}\notin{\mathcal{X}}.$ Therefore, with $x_{0}$ given
by (30), $x_{4}\notin{\mathcal{X}},$ and thus the assumptions of Theorem 4.5
are satisfied. In this case, Fig. 7 shows that the output does not converge
and is bounded, and thus the system has self-excited oscillations. In fact,
the asymptotic response is periodic with period 2 steps.

Figure 7: Example 5.1: a) and b) show the response of the Lur’e system in the
case where $x_{0}$ is given by (28) and c) and d) show the response of the
Lur’e system in the case where $x_{0}$ is given by (30). a) and c) show the
system response $y_{k}$ versus $k.$ b) and d) show the norm $||Px_{k}||$ of
the projection of $x_{k}$ onto ${\mathcal{R}}(\Psi)$ versus $k.$

##  6 Conclusions and Future Work

This paper analyzed a discrete-time Lur’e system that exhibits self-excited
oscillations. This system involves an asymptotically stable linear system with
a zero at 1 connected in feedback with a piecewise-linear saturation
nonlinearity. It was shown that, for sufficiently large loop gains, the
response of the system is bounded and does not converge, and thus the system
has self-excited oscillations. A numerical example illustrated the conditions
under which this discrete-time Lur’e system yields a self-excited response.
Future work will extend the Lur’e model to sigmoidal nonlinearities as well as
the use of this model for system identification.

## Appendix A

###### Lemma A.1

Let $z_{1},\ldots,z_{m}\in{\mathbb{C}}$ be distinct, assume that, for all
$i=1,\ldots,m,$ $|z_{i}|\geq 1,$ let $a_{1},\ldots,a_{m}\in{\mathbb{C}},$ and
assume that $\lim_{k\to\infty}\sum_{i=1}^{m}a_{i}z_{i}^{k}=0.$ Then,
$a_{1}=\cdots=a_{m}=0.$

Proof:  First, consider the case where $m=1.$ Suppose that $a_{1}\neq 0.$
Then, $\lim_{k\to\infty}a_{1}z_{1}^{k}=0,$ and thus
$\lim_{k\to\infty}|z_{1}|^{k}=0,$ which contradicts the assumption that
$|z_{1}|\geq 1.$ Therefore, $a_{1}=0.$

Next, consider the case where $m\geq 2.$ Let $s\in\\{1,\ldots,m\\},$ and
define the polynomial

$p_{s}(z)\stackrel{{\scriptstyle\triangle}}{{=}}\prod_{\genfrac{}{}{0.0pt}{}{1\leq
t\leq m}{t\neq s}}\frac{z-z_{t}}{z_{s}-z_{t}}.$

Writing $p_{s}(z)=\sum_{j=0}^{m-1}b_{j}z^{j},$ where
$b_{0},\ldots,b_{m-1}\in{\mathbb{C}}$, it follows that, for all $k\geq 1,$

$\sum_{i=1}^{m}a_{i}z_{i}^{k}p_{s}(z_{i})=\sum_{j=0}^{m-1}b_{j}\sum_{i=1}^{m}a_{i}z_{i}^{k+j}=\sum_{j=0}^{m-1}b_{j}w_{k+j},$
(32)

where
$w_{k}\stackrel{{\scriptstyle\triangle}}{{=}}\sum_{i=1}^{m}a_{i}z_{i}^{k}$.
Since, for all $j=0,\ldots,m-1,$ $\lim_{k\to\infty}w_{k+j}=0$, it follows that

$\lim_{k\to\infty}\sum_{i=1}^{m}a_{i}z_{i}^{k}p_{s}(z_{i})=0.$

Next, note that, since $p_{s}(z_{s})=1$ and, for all
$t\in\\{1,\ldots,m\\}\backslash\\{s\\},$ $p_{s}(z_{t})=0,$ it follows that

${\sum_{i=1}^{m}a_{i}z_{i}^{k}p_{s}(z_{i})}={a_{s}z_{s}^{k}}.$ (33)

It thus follows from (32) that

$a_{s}\lim_{k\to\infty}z_{s}^{k}=0.$ (34)

Now, suppose that $a_{s}\neq 0.$ Then, (34) implies that
$\lim_{k\to\infty}|z_{s}|^{k}=0.$ Hence, $|z_{s}|<1,$ which contradicts the
assumption that $|z_{s}|\geq 1.$ Therefore, $a_{s}=0,$ and thus
$a_{1}=\cdots=a_{m}=0.$ $\square$

The following is a corollary of Lemma A.1.

###### Corollary A.2

Let $z_{1},\ldots,z_{m}\in{\mathbb{C}}$ be distinct, assume that, for all
$i=1,\ldots,m,$ $|z_{i}|\geq 1,$ and let $a_{1},\ldots,a_{m}\in{\mathbb{C}}$,
at least one of which is nonzero. Then,

$\limsup_{k\to\infty}\left|\sum_{i=1}^{m}a_{i}z_{i}^{k}\right|>0.$ (35)

The following result is used in the proof of Lemma 3.2.

###### Proposition A.3

Let $\lambda_{1},\ldots,\lambda_{n}\in{\mathbb{C}}$ be distinct, and assume
that $\max_{i=1,\ldots,n}|\lambda_{i}|>1.$ Furthermore, let
$p_{1},\ldots,p_{n}$ be nonzero polynomials with complex coefficients, and,
for all $k\geq 1,$ define

$y_{k}\stackrel{{\scriptstyle\triangle}}{{=}}\sum_{i=1}^{n}p_{i}(k)\lambda_{i}^{k}.$
(36)

Then,

$\limsup_{k\to\infty}\left|y_{k}\right|=\infty.$ (37)

Proof:  Label $\lambda_{1},\ldots,\lambda_{n}$ such that
$|\lambda_{1}|\geq|\lambda_{2}|\geq\cdots\geq|\lambda_{n}|>0,$ and define
$\rho\stackrel{{\scriptstyle\triangle}}{{=}}\max_{i\in\\{1,\ldots,n\\}}\left|\lambda_{i}\right|>1$
and $s\in\\{1,\ldots,n-1\\}$ such that
$\rho=\left|\lambda_{1}\right|=\cdots=\left|\lambda_{s}\right|>\left|\lambda_{s+1}\right|.$
Furthermore, let $y_{k}=y_{k,1}+y_{k,2},$ where

$y_{k,1}\stackrel{{\scriptstyle\triangle}}{{=}}\sum_{i=1}^{s}p_{i}(k)\lambda_{i}^{k},\quad
y_{k,2}\stackrel{{\scriptstyle\triangle}}{{=}}\sum_{i=s+1}^{n}p_{i}(k)\lambda_{i}^{k},\\\
$ (38)

where, in the case $s=n,$ $y_{k,2}\stackrel{{\scriptstyle\triangle}}{{=}}0.$
Note that, for all $k\geq 0,$

$\displaystyle\left|\frac{y_{k,2}}{\rho^{k}}\right|$
$\displaystyle=\left|\sum_{i=s+1}^{n}p_{i}(k)\left(\frac{\lambda_{i}}{\rho}\right)^{k}\right|$
$\displaystyle\leq\left(\sum_{i=s+1}^{n}\left|p_{i}(k)\right|\right)\left|\frac{\lambda_{s+1}}{\rho}\right|^{k},$

which implies that

$\lim_{k\to\infty}\frac{y_{k,2}}{\rho^{k}}=0.$ (39)

Next, let $\lambda_{1},\ldots,\lambda_{s}$ be ordered such that
$d\stackrel{{\scriptstyle\triangle}}{{=}}\deg p_{1}=\deg p_{2}=\cdots=\deg
p_{m}>\deg p_{m+1}\geq\cdots\geq\deg p_{s},$ where $m\in\\{1,\ldots,s\\}.$
Thus, for all $i\in\\{1,\ldots,m\\},$ there exists a nonzero complex number
$a_{i}$ and a polynomial $q_{i}$ such that $\deg q_{i}<d$ and, for all $k\geq
0,$ $p_{i}(k)=a_{i}k^{d}+q_{i}(k).$ Hence,

$\frac{y_{k,1}}{\rho^{k}k^{d}}=\sum_{i=1}^{m}a_{i}z_{i}^{k}+\sum_{i=1}^{m}\frac{q_{i}(k)}{k^{d}}z_{i}^{k}+\sum_{i=m+1}^{s}\frac{p_{i}(k)}{k^{d}}z_{i}^{k},$
(40)

where, for all $i\in\\{1,\ldots,s\\},$
$z_{i}\stackrel{{\scriptstyle\triangle}}{{=}}\lambda_{i}/\rho.$ It thus
follows from (40) that

$\lim_{k\to\infty}\left(\frac{y_{k,1}}{\rho^{k}k^{d}}-\sum_{i=1}^{m}a_{i}z_{i}^{k}\right)=0.$
(41)

Hence, (39) and (41) imply that

$\lim_{k\to\infty}\left(\frac{y_{k}}{\rho^{k}k^{d}}-\sum_{i=1}^{m}a_{i}z_{i}^{k}\right)=0.$
(42)

Next, since $\lambda_{1},\ldots,\lambda_{m}$ are distinct, it follows that
$z_{1},\ldots,z_{m}$ are distinct. Corollary A.2 thus implies that there exist
$\varepsilon>0$ and an increasing sequence $(k_{j})_{j\geq 0}$ of positive
integers such that, for all $j\geq 0,$

$\left|\sum_{i=1}^{m}a_{i}z_{i}^{k_{j}}\right|>2\varepsilon,$ (43)

and, from (42), that there exists a nonnegative integer $j_{0}$ such that, for
all $j\geq j_{0},$

$\left|\frac{y_{k_{j}}}{\rho^{k_{j}}k_{j}^{d}}-\sum_{i=1}^{m}a_{i}z_{i}^{k_{j}}\right|<\varepsilon.$
(44)

Hence, (43) and (44) imply that, for all $j\geq j_{0},$

$\left|\frac{y_{k_{j}}}{\rho^{k_{j}}k_{j}^{d}}\right|\geq\left|\sum_{i=1}^{m}a_{i}z_{i}^{k_{j}}\right|-\left|\frac{y_{k_{j}}}{\rho^{k_{j}}k_{j}^{d}}-\sum_{i=1}^{m}a_{i}z_{i}^{k_{j}}\right|>\varepsilon,$

and thus $\left|y_{k_{j}}\right|>\varepsilon\rho^{k_{j}}k_{j}^{d}$. Since
$\rho>1$, it follows that

$\lim_{j\to\infty}\left|y_{k_{j}}\right|=\infty,$

which implies (37). $\square$

###### Lemma A.4

Let $M$ be an $n\times n$ real matrix, assume that 1 is not an eigenvalue of
$M,$ and let $x$ and $y$ be $n\times 1$ real vectors. Then the following
statements hold:

  1. $i)$

If $\lim_{k\to\infty}M^{k}$ exists, then the limit is zero.

  2. $ii)$

If $\lim_{k\to\infty}M^{k}x$ exists, then the limit is zero.

  3. $iii)$

If $\lim_{k\to\infty}y^{\rm T}M^{k}x$ exists, then the limit is zero.

Proof:  To prove iii), write $\chi_{M}(z)=z^{n}+a_{1}z^{n-1}+\cdots+a_{n}$,
and, for all $k\geq 0,$ define
$t_{k}\stackrel{{\scriptstyle\triangle}}{{=}}y^{\rm T}M^{k}x$. Then the
Cayley-Hamilton implies that, for all $k\geq 0,$

$t_{n+k}+a_{1}t_{n+k-1}+\cdots+a_{n}t_{k}=0.$ (45)

Since $\lim_{k\to\infty}t_{k}$ exists, letting $k\to\infty$ in (45) yields

$\displaystyle 0$ $\displaystyle=(1+a_{1}+\cdots+a_{n})\lim_{k\to\infty}t_{k}$
$\displaystyle=\chi_{M}(1)\lim_{k\to\infty}t_{k}.$ (46)

Since $\chi_{M}(1)\neq 0,$ (46) implies that $\lim_{k\to\infty}t_{k}=0.$
Finally, i) and ii) follow by taking $x$ and $y$ to be columns of the $n\times
n$ identity matrix. $\square$

## CRediT authorship contribution statement

Juan A. Paredes: Conceptualization, Formal analysis, Writing – original draft,
Writing – review & editing. Syed Aseem Ul Islam: Conceptualization, Writing –
original draft. Omran Kouba: Formal analysis, Writing – review & editing.
Dennis S. Bernstein: Writing – original draft, Writing – review & editing,
Supervision, Funding acquisition.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## Acknowledgements

The first author was supported by NSF grant CMMI 1634709.

## References

  * Ding [2010] W. Ding, Self-Excited Vibration: Theory, Paradigms, and Research Methods, Springer, 2010\.
  * Jenkins [2013] A. Jenkins,  Self-oscillation,  Physics Reports 525 (2013) 167–222.
  * Khalil [2002] H. K. Khalil, Nonlinear Systems, third ed., Prentice Hall, 2002\.
  * Mees and Chua [1979] A. Mees, L. Chua,  The hopf bifurcation theorem and its applications to nonlinear oscillations in circuits and systems,  IEEE Trans. Circ. Sys. 26 (1979) 235–254.
  * Tomberg and Yakubovich [1989] E. A. Tomberg, V. A. Yakubovich,  Conditions for auto-oscillations in nonlinear systems,  Siberian Mathematical Journal 30 (1989) 641–653.
  * Hang et al. [2002] C. Hang, K. Astrom, Q. Wang,  Relay feedback auto-tuning of process controllers—a tutorial review,  J. Process Control 12 (2002) 143–162.
  * Jian and Yu-shu [2004] X. Jian, C. Yu-shu,  Effects of time delayed velocity feedbacks on self-sustained oscillator with excitation,  Applied Math. Mech. 25 (2004) 499–512.
  * Stan and Sepulchre [2007] G. Stan, R. Sepulchre,  Analysis of interconnected oscillators by dissipativity theory,  IEEE Trans. Autom. Contr. 52 (2007) 256–270.
  * Aguilar et al. [2009] L. T. Aguilar, I. Boiko, L. Fridman, R. Iriarte,  Generating self-excited oscillations via two-relay controller,  IEEE Trans. Autom. Contr. 54 (2009) 416–420.
  * Chatterjee [2011] S. Chatterjee,  Self-excited oscillation under nonlinear feedback with time-delay,  J. Sound Vibr. 330 (2011) 1860–1876.
  * Risau-Gusman [2016] S. Risau-Gusman,  Effects of time-delayed feedback on the properties of self-sustained oscillators,  Phys. Rev. E 94 (2016) 042212.
  * Zanette [2017] D. H. Zanette,  Self-sustained oscillations with delayed velocity feedback,  Papers in Physics 9 (2017) 090003–1–090003–7.
  * Rasvan [1998] V. Rasvan,  Self-sustained oscillations in discrete-time nonlinear feedback systems,  in: Proc. 9th Mediterranean Electrotechnical Conf., 1998, pp. 563–565.
  * Gentile et al. [2011] F. S. Gentile, A. L. Bel, M. B. D’Amico, J. L. Moiola,  Effect of delayed feedback on the dynamics of a scalar map via a frequency-domain approach,  Chaos: An Interdisciplinary J. Nonlinear Science 21 (2011) 023117.
  * Nayfeh and Mook [2008] A. H. Nayfeh, D. T. Mook, Nonlinear oscillations, Wiley, 2008\.
  * Minorsky [2010] N. Minorsky,  Self-excited oscillations in dynamical systems possessing retarded actions,  in: R. Bellman, R. Kalaba (Eds.), Classic Papers in Control Theory, Dover, 2010, pp. 143–149.
  * Paredes et al. [2020] J. Paredes, S. A. Ul Islam, D. S. Bernstein,  A time-delayed Lur’e model with biased self-excited oscillations,  in: Proc. Amer. Contr. Conf., 2020, pp. 2699–2704.
  * Bernstein [2018] D. S. Bernstein, Scalar, Vector, and Matrix Mathematics: Theory, Facts, and Formulas—Revised and Expanded Edition, Princeton, 2018.
  * Meyer [2000] C. D. Meyer, Matrix Analysis and Applied Linear Algebra, SIAM, 2000.
  * Jiang and Wang [2001] Z.-P. Jiang, Y. Wang,  Input-to-state stability for discrete-time nonlinear systems,  Automatica 37 (2001) 857–869.