ArXiv preprint

[orcid=0000-0001-8013-8613<EMAIL_ADDRESS>]

<EMAIL_ADDRESS>]

<EMAIL_ADDRESS>]

<EMAIL_ADDRESS>]

# Fault Detection and Diagnosis with Imbalanced and Noisy Data: A Hybrid
Framework for Rotating Machinery

Masoud Jalayer  Amin Kaboli  Carlotta Orsenigo  Carlo Vercellis  Department of
Management, Economics and Industrial Engineering, Politecnico di Milano, Via
Lambruschini 24/b, 20156, Milan, Italy  Department of Mechanical Engineering,
University of Victoria, Victoria BC, V8P 5C2, Canada  Institute of Mechanical
Engineering, School of Engineering, Swiss Federal Institute of Technology at
Lausanne (EPFL), Lausanne, Switzerland

###### Abstract

Fault diagnosis plays an essential role in reducing the maintenance costs of
rotating machinery manufacturing systems. In many real applications of fault
detection and diagnosis, data tend to be imbalanced, meaning that the number
of samples for some fault classes is much less than the normal data samples.
At the same time, in an industrial condition, accelerometers encounter high
levels of disruptive signals and the collected samples turn out to be heavily
noisy. As a consequence, many traditional Fault Detection and Diagnosis (FDD)
frameworks get poor classification performances when dealing with real-world
circumstances. Three main solutions have been proposed in the literature to
cope with this problem: (1) the implementation of generative algorithms to
increase the amount of under-represented input samples, (2) the employment of
a classifier being powerful to learn from imbalanced and noisy data, (3) the
development of an efficient data pre-processing including feature extraction
and data augmentation. This paper proposes a hybrid framework which uses the
three aforementioned components to achieve an effective signal based FDD
system for imbalanced conditions. Specifically, it first extracts the fault
features, using Fourier and wavelet transforms to make full use of the
signals. Then, it employs Wasserstein Generative Adversarial Networks (WGAN)
to generate synthetic samples to populate the rare fault class and enhance the
training set. Moreover, to achieve a higher performance a novel combination of
Convolutional Long Short-term Memory (CLSTM) and Weighted Extreme Learning
Machine (WELM) is also proposed. To verify the effectiveness of the developed
framework, different bearing datasets settings on different imbalance
severities and noise degrees were used. The comparative results demonstrate
that in different scenarios GAN-CLSTM-ELM significantly outperforms the other
state-of-the-art FDD frameworks.

###### keywords:

Surface Inspection Optical Quality ControlComputer VisionImage
AugmentationImage Object DetectionFault Diagnosis

##  1 Introduction

Rotating machinery is one of the essential equipment in today’s industrial
environments. From petroleum, automobile, chemicals, pharmaceutical, mining,
power generation plants to consumer goods, at least there is a machine with a
rotating component. The rotating component could be the gearbox, axles, wind,
steam and gas turbines, centrifugal and oil-free screw compressors, and pumps.
30% of rotating machinery breakdowns are mainly caused by loose, partially
rubbed, misaligned, cracked, and unbalanced rotating parts [1]. Machine
breakdowns can present complex challenges during day-to-day operations and
significantly impact business profitability and operations productivity.
Monitoring machine health conditions can prevent machine breakdowns and reduce
the maintenance costs of manufacturing systems [2]. It is, hence, crucial to
develop efficient diagnosis systems to analyze different health conditions of
the rotating components.

Figure 1: Main steps of an automated Fault Detection and Diagnosis system

There are two main approaches for coping with fault detection and diagnosis in
rotating machinery: (1) physical-based control systems and (2) data-driven-
based models. Recent advancements in computer processing and digital
technologies enhanced the robustness and higher computational capabilities to
use data-driven fault detection and diagnosis models. Implementing these
models enable us to monitor and control the parameters of machines from a
remote distance and drive insights. That is the main reason for which data-
driven fault detection and diagnosis models are used in smart manufacturing
systems [2].

The main contributions of this paper are as follows: (1) In order to get
higher classification performance in different environments, a hybrid deep
learning architecture is designed such that it takes Fourier and Wavelet
spectra of the vibration signals. This architecture uses CNN blocks to find
shift-agnostic characteristics of the fault types, a LSTM block which
understands the spatiotemporal and sequential features of it and, finally, a
Weighted ELM classifier which is effective in learning from scarce patterns,
the necessity of which is examined through experimental comparisons. The
proposed classifier is named CLSTM-ELM. (2) A Wasserstein-GAN model with a
gradient penalty is developed and employed in the hybrid framework to
reproduce rare patterns and enhance the training set. The effectiveness of
this proposition is investigated in Section 5. (3) A comprehensive set of
scenarios is designed to study the effect of different imbalance severities
and noise degrees on the performance of the framework. A sensitivity analysis
is conducted on the scenarios revealing more insights about the
characteristics of the model. (4) Seven state-of-the-art FDD models are chosen
to compete with the proposed framework on four different dataset settings. The
experimental comparison illustrates how implementing WGAN-GP and W-ELM can
improve the classifier performance and shows the superiority of GAN-CLSTM-ELM
over other algorithms.

The rest of the paper is organized as follows. Section 2 provides an overview
of the principal AI-based approaches proposed for FDD problems. In Section 3,
the theory behind WGAN-GP, LSTM, Convolutional layers and W-ELM is briefly
reviewed. Then, the proposed hybrid framework, GAN-CLSTM-ELM, is presented in
Section 4. Section 5 compares the performance of different FDD algorithms on
different imbalance ratios and noise severities. Finally, some research
conclusions and future extensions are provided in Section 6.

##  2 Review of Current Models

Early data-driven fault detection and diagnosis (hereafter FDD) models have
benefited from traditional Artificial Intelligence (AI) models, or “shallow
learning” models, such as Support Vector Machines (SVM), Decision Trees (DT),
and Multi-layer Perceptron (MLP) [3]. Despite the applicability of traditional
AI models to FDD problems, these models show poor performances and limitations
when dealing with complicated fault patterns such as the above-mentioned
rotating machinery faults [4]. One of the first applications of rotating
machinery FDD dates back to 1969 in Boeing Co., when Balderston [5]
illustrated some characteristics of the fault signs on the signals measured by
an accelerometer in natural and high frequencies. [6] employed the rectified
envelope signals with a synchronous averaging, which was later called
“envelope analysis”, to identify bearing local faults. The peak localization
in the vibration signal spectrum is another classical example of the fault
detection methods for the ball bearing faults [7].

Recently, with the emergence of novel deep learning (DL) architectures and
their promising pattern recognition capabilities, many researchers proposed
deep learning solutions for data-driven-based FDD systems [8]. These FDD
approaches rely on the common assumption that the distribution of classes for
different machine health conditions is approximately balanced. In practice,
however, the number of instances may significantly differ from a fault class
to another. This causes a crucial issue since a classifier which has been
trained on such a data distribution primarily exhibits a skewed accuracy
towards the majority class, or fails to learn the rare patterns. Most of the
proposed FDD approaches, thus, suffer from higher misclassification ratios
when dealing with scarce conditions such as in high-precision industries where
the number of faults are limited [9].

Through their deep architectures, DL-based methods are capable of adaptively
capturing the information from sensory signals through non-linear
transformations and approximate complex non-linear functions with small errors
[3]. Auto-encoders (AE) are among the most promising DL techniques for
automatic feature extraction of mechanical signals. They have been adopted in
a variety of FDD problems in the semiconductor industry [10], foundry
processes [11], gearboxes [12] and rotating machinery [13], [14]. [15]
employed the “stacked” variation of AE to initialize the weights and offsets
of a multi-layer neural network and to provide an expert knowledge for
spacecraft conditions. However, to cope with mechanical signals, using a
single AE architecture has shown some drawbacks: it may only learn similar
features in feature extraction and the learned features may have shift variant
properties which potentially lead to misclassification. Some approaches were
proposed to make this architecture appropriate for signal-based fault
diagnosis tasks. [16] used a local connection network on a normalized sparse
AE, called NSAE-LCN, to overcome these shortcomings. [17] developed a stacked-
AE to directly learn features of mechanical vibration signals on a motor
bearing dataset and a locomotive bearing dataset; specifically, they first
used a two-layer AE for sparse filtering and then applied a softmax regression
to classify the motor condition. The combination of these two techniques let
the method achieved high accuracy in bearing fault diagnosis.

Extreme learning machine (ELM) is a competitive machine learning technique,
which is simple in theory and fast in implementation. As an effective and
efficient machine learning technique, ELM has attracted tremendous attention
from various fields in recent years. Some researchers suggest ELM and Online
Sequential ELM (OS-ELM) for learning from imbalance data [18]; [19]; [20]. ELM
and OS-ELM can learn extremely fast due to their ability to learn data one-by-
one or chunk-by-chunk [21]. Despite their effective performances on online
sequential data, the performance associated to their classical implementation
on highly imbalanced data is controversial; according to [22], for example,
OS-ELM tends to have poor accuracy on such data. Therefore, they proposed a
voting-based weighted version of it, called VWOS-ELM, to cope with severely
rare patterns, whereas [9] developed a two-stage hybrid strategy using a
modified version of OS-ELM, named PL-OSELM. In offline stage, the principal
curve is employed to explore the data distribution and develop an initial
model on it. In online stage, some virtual samples are generated according to
the principal curve. The algorithm chooses virtual minority class samples to
feed more valuable training samples.

Considering the promising results obtained by ELM- based classifiers coping
with imbalanced data, they accordingly became one of the mainstreams in FDD
research area. In [23], the authors developed an evolutionary OS-ELM for FDD
for bearing elements of high-speed electric multiple units. They employed a
K-means synthetic minority oversampling technique (SMOTE) for oversampling the
minority class samples. They also used an artificial bee colony (ABC)
algorithm to find a near-optimum combination of input weights, hidden layer
bias, and number of hidden layer nodes of the OS-ELM. In another paper, [24]
used density-weighted one-class ELM for fault diagnosis in high-voltage
circuit breakers (HVCBs), using vibration signals. [25] applied an adaptive
class-specific cost regulation ELM (ACCR-ELM) with variable-length brainstorm
algorithm for its parameter optimization to conveyor belt FDD. The proposed
algorithm exhibits a stable performance under different imbalance ratios. [26]
presented a feature extraction scheme on time-domain, frequency-domain, and
time-frequency-domain, to feed a full spectrum of information gained from the
vibration signals to the classifier. They also demonstrated that the cost-
sensitive gradient boosting decision tree (CS-GBDT) shows a satisfactory
performance for imbalanced fault diagnosis. In another FDD framework for
rolling bearings [27], the authors coupled an Optimized Unsupervised Extreme
Learning Machine (OUSELM) with an Adaptive Sparse Contractive Auto-encoder
(ASCAE). The ASCAE can gain an effectual sparse and more sensitive feature
extraction from the bearing vibration signals. A Cuckoo search algorithm was
also proposed to optimize the ELM hyper-parameters. Another variation of ELM
was developed by [28] to deal with imbalanced aircraft engines fault data
which is derived from the engine’s thermodynamic maps. This ELM variation
flexibly sets a soft target margin for each training sample; hence, it does
not need to force the margins of all the training samples exactly equaling one
from the perspective of margin learning theory. After some experiments on
different datasets, including the aircraft engine, it is concluded that SELM
outperforms ELM.

On the other hand, there are frameworks for imbalanced and noisy FDD without
the employment of any ELM variation. [16] proposed a Deep Normalized
Convolutional Neural Network (DNCNN) for FDD under imbalanced conditions. The
DNCNN employs a weighted softmax loss which assumes that the misclassification
errors of different health conditions share an equivalent importance.
Subsequently, it minimizes the overall classification errors during the
training processes and achieves a better performance when dealing with
imbalanced fault classification of machinery adaptively. [29] used WGAN-GP to
interpolate stochastically between the actual and virtual instances so that it
ensures that the transition region between them is stable. They also utilized
a Stacked-AE to classify the enhanced dataset and determined the availability
of the virtual instances. Since a single GAN model encounters hardship and
poor performance when dealing with FDD datasets, [30] proposed a framework
based on GANs under small sample size conditions which boost the adaptability
of feature extraction and consequently diagnosis accuracy. The effectiveness
and satisfactory performance of the proposed method were demonstrated using
CWRU bearing and gearbox datasets. Another novel GANs-based framework, named
dual discriminator conditional GANs (D2CGANs), has been recently proposed to
learn from the signals on multi-modal fault samples [31]. This framework
automatically synthesizes realistic high-quality fake signals for each fault
class and is used for data augmentation such that it solves the imbalanced
dataset problem. After some experiments on the CWRU bearing dataset, the
authors showed that Conditional-GANs, Auxiliary Classifier-GANs and D2CGANs
significantly outperform GANs and Dual-GANs. [32] proposed a framework which
adopts a CNN-based GANs with the coordinative usage of two auxiliary
classifiers. The experimental results on analog-circuit fault diagnosis data
suggested that the proposed framework achieves a better classification
performance than that of DBN, SVM and artificial neural networks (ANN). [33]
presented a CNN-based GANs for rotating machinery FDD which uses a Wavelet
Transform (WT) technique. The proposed so-called WT-GAN-CNN approach extracts
time-frequency image features from one-dimension raw signals using WT.
Secondly, GANs are used to generate more training image samples while the
built CNN model is used to accomplish the FDD on the augmented dataset. The
experimental results demonstrated high testing accuracy in the interference of
severe environment noise or when working conditions were changed.

##  3 Background Theory

###  3.1 WGAN-GP

In the classical definition, GANs consist of two adversarial networks trained
in opposition to one another: (1) a generative model, $G$, which learns the
data distribution to generate a fake sample $\tilde{x}^{(i)}=G(z)$ from a
random vector $z$, where $z\sim\mathscr{P}_{z}$ and $\mathscr{P}_{z}$ is the
noise distribution; (2) a discriminator model, $D$, which determines if the
sample is generated by $G$ or is a real sample. $G$ strives to deceive $D$ by
making realistic random samples, while $D$ receives both real and fake
samples. On the contrary, $D$ tries to find out the source of each sample by
calculating its corresponding probability, $p(S|x)=D(x)$, and is trained to
maximize the log-likelihood it assigns to the correct source [34]

$\begin{split}\min_{G}\max_{D}V(D,G)=\mathbb{E}_{x\sim\mathscr{P}_{r}}[\log(D(x))]+\mathbb{E}_{\tilde{x}\sim\mathscr{P}_{f}}[\log(1-D(\tilde{x}))],\end{split}$
(1)

where $\mathscr{P}_{r}$ and $\mathscr{P}_{f}$ denote the distribution of the
raw data and of the fake samples, respectively. Then, the model reaches a
dynamic equilibrium if $\mathscr{P}_{f}=\mathscr{P}_{r}$ [34].

While GAN is a powerful generative model, it suffers from training
instability. Some different solutions have been proposed to solve this
problem. Wasserstein GAN (WGAN) is one of the novel proposed techniques
offering a new loss function which has demonstrated a better performance and a
better model stability. The present paper uses a variation of GAN, Entropy-
based WGAN-GP proposed by [35], generating an entropy-weighted label vector
for each class with respect to its frequency.

When the discriminator $D$ is sufficiently trained, the gradient of the
generator $G$ is relatively small; when the effect of $D$ is lower, it gets
larger. WGAN employs a distance called Wasserstein to calculate the difference
between the distributions of the real and fake samples [36], which can be
mathematically written as follows:

$\mathcal{W}(\mathscr{P}_{r},\mathscr{P}_{f})=\inf_{\lambda\in\Pi(\mathscr{P}_{r},\mathscr{P}_{f})}\mathbb{E}_{(a,b)\sim\lambda}[\|a-b\|],$
(2)

where $\Pi(\mathscr{P}_{r},\mathscr{P}_{f})$ denotes the set of all
distributions with margins of $\mathscr{P}_{r}$ and $\mathscr{P}_{f}$, and
$\lambda(a,b)$ represents the distance between two given distributions, $a$
and $b$. Therefore, the Wasserstein variable can be interpreted as the
transportation cost between the distributions of real and fake datasets. To
avoid the gradient uninformativeness issue and to guarantee the existence and
uniqueness of the optimal discriminative function and the respective Nash
equilibrium, Lipschitz condition is applied [37]. In the proposed WGAN, the
discriminative function is restricted to 1-Lipschitz. The WGAN, hence,
proposes a revised objective function based on 1, using Kantorovich-Rubinstein
duality, and is formulated as:

$\min_{G}\max_{D\in\mathcal{L}_{1}}\mathbb{E}_{x\sim\mathscr{P}_{f}}[D(x)]-\mathbb{E}_{\tilde{x}\sim\mathscr{P}_{r}}[D(\tilde{x})],$
(3)

where $\mathcal{L}_{1}$ is the collection of 1-Lipschitz functions. In this
case, under an optimal discriminator which minimizes the objective function
with respect to the parameters of $G$, the model strives to minimize
$\mathcal{W}(\mathscr{P}_{r},\mathscr{P}_{f})$.

Let $\delta\sim U[0,1]$ be a random number to have a linear interpolation
between $x$ and $\tilde{x}$:

$\hat{x}=\delta x+(1-\delta)\tilde{x},$ (4)

Therefore, the loss function of the discriminator, $\mathscr{L}_{D}$, can be
determined as follows

$\mathscr{L}_{D}=\mathbb{E}_{\tilde{x}\sim\mathscr{P}_{g}}[D(\tilde{x})]-\mathbb{E}_{x\sim\mathscr{P}_{r}}[D(x)]+\gamma\mathbb{E}_{\hat{x}\sim\mathscr{P}_{z}}[(\|\nabla_{\hat{x}}D(\hat{x})\|_{2}-1)^{2}],$
(5)

where $\gamma$ stands for the gradient penalty coefficient. The last part of
Eq.5 denotes the gradient penalty:
$\gamma\mathbb{E}_{\hat{x}\sim\mathscr{P}_{z}}[(\|\nabla_{\hat{x}}D(\hat{x})\|_{2}-1)^{2}]$
[38].

Figure 2 exhibits the simplified process of synthesizing rotating machinery
signal samples out of some random noises in a conventional GAN model.

Figure 2: The schematic process in GANs

###  3.2 CLSTM

RNN (Recurrent Neural Network) is a powerful class of artificial deep learning
architectures proposed to identify patterns in sequential data. It can
consider time and sequence by taking the present data and the recent past data
as inputs. RNN is trained across the time steps using backpropagation.
However, due to the multiplication of gradients at time steps the gradient
value may vanish or blow up rapidly. This issue limits its usage when the time
window is greater than 10 discrete time steps [39]. By adding constant error
carousels and introducing forget gates to RNN, a new form of RNN architecture
is proposed, named LSTM. These adopted forget gates are able to control the
utilization of information in the cell states and impede the vanishing or
exploding gradient issues [40]. Compared to RNN, LSTM is more powerful in
capturing long-term dependencies of the data features where it can handle
time-windows exceeding 1000 discrete time stamps [39].

On the other hand, convolutional neural networks are mainly composed of
convolutional, pooling and normalization layers, making them capable of
understanding the data shift-invariance and sharing the weights through
convolutional connections. This weight sharing makes CNN lighter for
computation since it reduces the number of parameters in the network.

Let $x_{i}=[\kappa_{1},\ldots,\kappa_{L}]$ be the sequential data, where $L$
denotes the length of the signal sample, and $\kappa_{i}\in\mathbb{R}^{d}$ the
set of values at each timestamp, where $d$ is the number of channels. The
convolution operation is defined as the following dot product:

$c_{i}=\varphi(u\cdot\kappa_{i:i+m-1}+b),$ (6)

where $u\in\mathbb{R}^{md}$ is a filter vector, $b$ is the bias, $\varphi$ is
an activation function, and $\kappa_{i:i+m-1}$ represents an $m$-length window
starting from the $i_{th}$ time stamp of the sample. Sliding the filter window
from the first timestamp of the sample to its last possible one, a feature map
is given as follows:

$\mathscr{C}_{\xi}=[c_{1},c_{2},\dotsc,c_{L-m+1}],$ (7)

where $\mathscr{C}_{\xi}$ corresponds to the $\xi_{th}$ filter.

To reduce the length of these feature maps and minimizing the model
parameters, pooling layers are proposed. Max Pooling layer is one of the most
common pooling techniques. The compressed feature vector, $\mathscr{H}$, is
defined as follows:

$\mathscr{H}=[h_{1},h_{2},\dotsc,h_{(\frac{L-m}{s})+1}],$ (8)

where
$h_{\xi}=\max(\mathscr{C}_{(\xi-1)s},\mathscr{C}_{(\xi-1)s+1},\dotsc,\mathscr{C}_{(\xi
s-1)})$ on the $s$ consecutive values of feature map $\mathscr{C}_{\xi}$ and
$s$ denotes the pooling length. Batch normalization is widely used in CNN
blocks to reduce the shift of interval covariance and make the learning
quicker by alleviating the computational load. The normalization is completed
by making each individual scalar feature with zero mean and unit variance.
This process can be mathematically described as:

$\hat{\kappa_{i}}=(\kappa_{i}-\mathbb{E}[x_{i}])\sqrt{Var(x_{i})+\epsilon},$
(9)

where $\epsilon$ is a small constant added for numerical stability. However,
the extracted features can be affected when the features of a certain layer
are normalized directly by Eq. 9, leading to poor network expression
abilities. To resolve this issue, each normalized value $\kappa_{i}$ is
modified based on the scale parameter $\varrho_{i}$ and the shift parameter
$\varpi_{i}$. These two learnable reconstruction parameters can recover the
feature distribution of the original network. The following formula can be
used to determine the output of the neuron response:

$\hat{\nu_{i}}=\varrho_{i}\hat{\kappa_{i}}+\varpi_{i}.$ (10)

###  3.3 W-ELM

ELM is in general a single-hidden-layer feed-forward neural network (SLFN).
The difference between SLFN and ELM lies in how the weights of hidden layer
and output layer neurons are updated. In SLFN, the weights of both input and
outputs layers are initialized randomly, and the weights of both the layers
are updated by the backpropagation algorithm. In ELM, the weights of the
hidden layers are assigned randomly but never updated, and only the weights of
the output layer are updated during the training process. As in ELM, the
weights of only one layer are to be updated as opposed to both layers of SLFN;
this makes ELM faster than SLFN.

Let the training dataset be $\\{(x_{i},y_{i})\\}_{(i=1)}^{N}$, where $x_{i}$
is the input vector and $y_{i}$ is the output vector. The output of the
$j^{th}$ hidden layer neuron is given by $\varphi_{a}(a_{j},b_{j},x_{i})$,
where $a_{j}$ is the weight vector connecting the input neurons to the
$j^{th}$ hidden layer neuron, $b_{j}$ is the bias of the $j^{th}$ hidden
neuron, and $\varphi_{a}$ is the activation function. Each hidden layer neuron
of ELM is also connected to each output layer neuron with some associated
weights. Let $\beta=[\beta_{1},\beta_{2},\dotsc,\beta_{K}]^{T}$ denote the
output weights connecting the hidden layer (composed of $K$ hidden nodes) with
output neurons. Thus, the $i^{th}$ output is determined as:

$o_{i}=\sum^{K}_{j=1}\beta_{j}\phi_{a}(a_{j},b_{j},x_{i}),\quad i=1,\dots,N.$
(11)

Let $H=(H_{ij})=(\phi_{a}(w_{j},b_{j},x_{i}))$ be the hidden layer matrix. The
$N$ equations of the output layer (Eq. 11) can be shortly written as follows:

$O=H\beta.$ (12)

Using Moore-Penrose generalized inverse [41], $H^{\dagger}$, a least square
solution, referred to as the extreme learning machine, can be determined
mathematically as follows:

$\beta=H^{\dagger}Y=\begin{cases}H^{T}(\frac{1}{C}+HH^{T})^{-1}Y&N<K\\\
(\frac{1}{C}+H^{T}H)^{-1}H^{T}Y&N\geq K,\end{cases}$ (13)

where $C$ is a positive parameter to achieve a better generalization
performance [42]. Weighted ELM [20] considers a $N$×$N$ diagonal matrix $W$
associated with each training sample $x_{i}$. If $x_{i}$ belongs to the
minority class, the algorithm allocates a relatively larger weight to $w_{i}$
rather than those of majority classes, which intensifies the impact of
minority classes in the training phase. Therefore, the solution of Eq. 12 will
be obtained by using the optimization formula of ELM:

$\begin{split}\emph{Minimize}:L_{P^{ELM}}=\frac{1}{2}\|\beta\|^{2}+CW\frac{1}{2}\sum^{N}_{i=1}\|\eta_{i}\|^{2}\\\
\emph{subject
to}:h(x_{i})\beta=y^{T}_{i}-\eta^{T}_{i},i=1,\dotsc,N.\end{split}$ (14)

According to the KKT theorem [43], the solution to Eq. 14 is obtained as
follows:

$\beta=H^{\dagger}Y=\begin{cases}H^{T}(\frac{1}{C}+WHH^{T})^{-1}WY&N<K\\\
(\frac{1}{C}+H^{T}WH)^{-1}WH^{T}Y&N\geq K.\end{cases}$ (15)

##  4 The Proposed FDD Model

As mentioned in Section 1, there is a need to improve the performance of
imbalanced and noisy FDD systems. Therefore, this paper presents a hybrid
framework which embeds three steps: (1) the employment of a generative
algorithm to improve the training set, (2) the signal processing with FFT and
CWT techniques providing the DL classifier a deeper understanding of the
fault’s identity, (3) the development of a hybrid classifier based on CNN,
LSTM and weighted ELM, as illustrated in Figure 3.

###  4.1 Sample generation model design

The structure of the WGAN-GP generator $G$ comprises a five-layered
autoencoder of $l,\frac{l}{2},\frac{l}{4},\frac{l}{2}$ and $l$ neurons, while
the discriminator $D$ is composed of three convolutional-LSTM blocks. The
input variable $z$ has a dimension of $l\times 1$. Due to the poor performance
in weight clipping of WGAN-GP, the paper uses an alternative in the form of a
gradient penalty in the discriminator loss function, which has been introduced
by [44] and which achieves high performances compared to other GAN models.
Algorithm 1 shows how WGAN-GP works, where $\rho$ is the gradient coefficient,
$\chi_{d}$ is the number of discriminator iterations per each generator
iteration, $\chi_{b}$ denotes the batch size, $\vartheta$,$\mu_{1}$ and
$\mu_{2}$ are the Adam hyperparameters, and $\omega_{1}$ and $\theta_{1}$
represent the initial discriminator and generator parameters, respectively.

Input: $\gamma,n_{critic},m,\alpha,\beta_{1},\beta_{2},\omega_{0},\theta_{0}$

while  _$\theta$ is not converged_ do

for  _$t\leftarrow 1,...,n_{critic}$_ do

for  _$i\leftarrow 1,...,m$_ do

Sample from real dataset $x\sim\mathscr{P}_{r}$,

Generate noise samples $z\sim\mathscr{P}_{z}$,

Generate a random number $\delta\sim U[0,1]$

$\tilde{x}\leftarrow G_{\theta}(z)$

$\hat{x}\leftarrow\delta x+(1-\delta)\tilde{x}$

$\mathscr{L}^{(i)}\leftarrow
D_{\omega}(\tilde{x})-D_{\omega}(x)+\gamma(\|\nabla_{\tilde{x}}D_{\omega}(\hat{x})\|_{2}-1)^{2}$

$\omega\leftarrow
Adam(\nabla_{\omega}\frac{1}{m}\sum_{i=1}^{m}L^{i},\omega,\alpha,\beta_{1},\beta_{2})$

Sample batch of $m$ noise samples $\\{z^{i}\\}^{m}_{i=1}\sim\mathscr{P}_{z}$

$\theta\leftarrow
Adam(\nabla_{\theta}\frac{1}{m}\sum_{i=1}^{m}-D_{\omega}(G_{\theta}(z)),\omega,\alpha,\beta_{1},\beta_{2})$

Algorithm 1 WGAN-GP

###  4.2 Fault diagnosis model design

In order to reveal more information about the fault characteristics, the paper
separately employs two signal processing feature extraction techniques on the
input samples: FFT and CWT [45], whose results are merged to be passed to the
deep learning architectures. As it is demonstrated in [46], employing these
two feature extraction techniques significantly improves the diagnosis
performance and increases the accuracy.

As it is shown in Figure 3, the paper proposes a dual-path deep learning
architecture which combines LSTM and CNN. The reason behind this duality lies
in the nature of these two architectures and the fact that each of them
explains a different feature type. More specifically, [47] illustrated that
the concatenation of CNN and LSTM features meaningfully enhances the
classification accuracy.

In the first pathway, after applying a one-dimensional convolutional layer on
the pre-processed input tensors which extracts the local and discriminative
features, an LSTM is added to encode the long-term temporal patterns. The
importance of adding a convolutional layer prior to the LSTM layer is not only
that it reduces the high-frequency noise impact, but it also helps the LSTM to
learn more effectively. The convolutional layer processes the sequence of the
features extracted after FFT and CWT. The model slides the kernel filters on
the sequence and generates feature maps. These feature maps are subsequently
processed by an LSTM which acquires the spatial dependencies of these
sequenced features.

On the other hand, in the second pathway three one-dimensional CNN blocks are
adjoined to better extract the local features of the Fourier and Wavelet
transform-based diagrams. Each CNN block contains a convolutional layer which
convolves the local regions of the diagrams, and a Rectified Linear Unit
(ReLU) activation function, which helps the network achieve a non-linear
expression and, consequently, make the learned features more distinguishable.
To reduce the computational complexity and decrease the covariance of shift
intervals, a batch normalization layer is added to each CNN block, following
by a max pooling layer which compresses the learnt features, advances its
local translation invariance and also alleviates the learning computational
expenses. These CNN blocks are followed by a flatten layer that reshapes the
tensor size to a vector to become compatible to join the output of the first
pathway and to be fed to the classifier.

For the classification architecture, the paper employs a weighted ELM which is
introduced in Section 3.3. ELM classifiers can get fast training speeds by
means of non-tuned training strategy. They also tend to have high
generalization performance in multi-class problems and showed excellent
classification performance in different studies [19]. Compared to unweighted
ELM, weighted ELM is aimed to put an additional accent on the samples implying
the imbalanced class distribution, so that the features in samples from the
minority class are also well perceived by the ELM [20]. Therefore, after
concatenating the outputs of both pathways, their combined learnt features are
passed to the W-ELM to diagnose the fault types.

###  4.3 General procedure of the proposed model

The schematic flowchart of the proposed intelligent FDD model is illustrated
in Figure 3. The general procedure of this model is summarized as follows:

  * •

Step 1: The sensory signals are collected from the accelerometers mounted on
the rotating machinery.

  * •

Step 2: The training, the test, and the validation datasets are constructed
from the raw signals to separate bursts by resampling.

  * •

Step 3: The training dataset is augmented using WGAN-GP introduced in Section
3.1 on the minority classes. The fake samples are added to the real samples to
make the training dataset balanced.

  * •

Step 4: By employing FFT and CWT techniques the model can extract fault
signatures which were hidden in the raw signals. The extracted Fourier and
Wavelet transform-based diagrams are concatenated to form three-dimensional
matrices which can be given in input to the deep learning blocks.

  * •

Step 5: These pre-processed samples go through two different paths of deep
learning blocks: (1) a one -dimensional convolutional layer followed by an
LSTM block, and (2) three blocks of CNN architectures followed by flatten and
dense layers.

  * •

Step 6: After concatenating the outputs of the two deep learning paths, a
W-ELM technique is used to classify the extracted deep features and diagnose
the fault type.

Figure 3: Schematic illustration of the proposed model for the training set

##  5 Results

To evaluate the effectiveness of the proposed method, some experiments were
run on one of the most widely used bearing fault datasets, known as Case
Western Reserve University (CWRU) bearing
dataset111https://csegroups.case.edu/bearingdatacenter/home. To conduct a
comprehensive comparison, we defined different noise and imbalance conditions
on which eight different DL-based FDD methods were tested. All experiments
were performed by using Python 3.9 on a computer with a GPU of NVIDIA Geforce
GTX 1070 with CUDA version of 10.1 and 16 GB of memory.

###  5.1 Dataset description

The paper employs CWRU bearing dataset using the test stand shown in Figure 4,
that consists of a motor, a torque transducer/encoder, a dynamometer, and
control electronics. The dataset consists of five different types of faults
corresponding to inner race, balls and outer race in three different
orientations: 3 o’clock (directly in the load zone), 6 o’clock (orthogonal to
the load zone) and 12 o’clock (opposite to the load zone). Moreover, the
faults are collected in a range of severity varying between 0.007 inches to
0.040 inches in diameter. The dataset is also recorded for motor loads, from 0
to 3 horsepower. However, for the sake of simplicity this paper uses only one
motor speed of 1797 RPM. The samples are collected at 12,000 samples/second
frequency from two accelerometers mounted on fan-end and drive-end of the
machine. In the experiments we took signal bursts of 800 timestamps, equal to
66.6 milli-seconds, to generate some different datasets of approximately
25,500 signal bursts.

To explore the diagnostic capabilities of the proposed framework in imbalanced
conditions, some non-equitant sets of samples were selected such that a fault
class becomes rare. Table 1 shows the distribution of samples for each machine
condition in the selected sets, where the value of $\alpha$ denotes the
percentage of minority class within the whole dataset. Accordingly, as
$\alpha$ decreases the imbalance degree increases. In this paper we chose
“out3” class to represent the minority class, whose samples correspond to the
outer race faults of opposite load zone position. In these scenarios, the
“health” class, corresponding to the healthy condition, represents
($80-\alpha$) percent of the whole dataset, while the other fault classes
account for 5% each. The generative algorithm, subsequently, strives to
equalize the sample size of the fault classes in the training set by
augmenting the minority class.

Adding “additive white Gaussian noise” with different signal-to-noise ratios
(SNRs) to the original samples, the paper is able to examine the performance
of GAN-CLSTM-ELM framework on different natural noise severities. These noisy
samples better portray the real-world industrial production settings where the
noise varies a lot. The original drive-end and fan-end signals with their
driven noisy samples are exhibited in Figure 5.

Figure 4: Two-horsepower (left), a torque transducer and encoder (center) and a dynamometer (right) used to collect the dataset Figure 5: Some noisy signal samples generated from raw sensory data with different SNRs Table 1: The distribution of condition type samples in the cases _minority share(%)_ | _Percentage of training samples in each condition_
---|---
_health_ | _inner_ | _ball_ | _out1_ | _out2_ | _out3_
$\alpha=4$ | 76% | 5% | 5% | 5% | 5% | 4%
$\alpha=2$ | 78% | 5% | 5% | 5% | 5% | 2%
$\alpha=1$ | 79% | 5% | 5% | 5% | 5% | 1%
$\alpha=0.5$ | 79.5% | 5% | 5% | 5% | 5% | 0.5%
$\alpha=0.25$ | 79.75% | 5% | 5% | 5% | 5% | 0.25%

###  5.2 GAN model selection

Figure 6: The generator and discriminator loss values for different GAN
architectures

As it is mentioned in the previous section, the proposition of Wasserstein
loss function and adding the gradient penalty to its loss function help
stabilize the generative algorithm. Figure 6 depicts how the proposed WGAN-GP
reaches an equilibrium after 9000 epochs where it can generate realistic
samples. Whereas, the other GAN generators make samples which cannot devise
their discriminators. As it can be clearly seen in Figure 6 their generator
loss values go significantly higher than those of the discriminators. This
comparison demonstrates why the implementation of WGAN-GP is preferred. Figure
7 shows some real samples of normal baseline, and fault conditions associated
with the bearing ball, inner race and outer race with fault diameters of 7
mils and 21 mils. Figure 8, similarly, visualizes the synthetic samples
generated by WGAN-GP after 10000 epochs.

Figure 7: Some real samples associated with different running conditions
Figure 8: Some random synthesized samples associated with different running
conditions made by WGAN-GP

###  5.3 The sensitivity analysis

In this section the paper illustrates a sensitivity analysis on the
performance of the proposed model by changing the $\alpha$ values and the
SNRs. Specifically, we considered 25 points with respect to
$\alpha={2^{k}:k=-2,-1,0,1,2}$ and SNR $=(10,20,50,75,100)$, and run the model
10 times at each point to achieve a robust analysis. Figure 9 and Figure 10
demonstrate the performance of GAN-CLSTM-ELM model with different metrics on
these points.

Figure 9: Accuracy and Recall performances of GAN-CLSTM-ELM in different SNR
and $\alpha$ levels Figure 10: AUC and $f_{1}$ score performances of GAN-
CLSTM-ELM in different SNR and $\alpha$ levels

As it can be seen in the figures, high levels of noise impact on the
performance of the model, changing the $f_{1}$ score from 100% to 95.91%, and
from 99.7% to 81.45% when $\alpha=4$ and $\alpha=0.25$, respectively. In this
defined space the accuracy, AUC and recall values fall above 96.7%, 92.6% and
81.16%, respectively. The model shows a relatively high robustness to both
noise and imbalance severities for SNRs greater than 20. At its best-case
scenario, where $\alpha=4$ and SNR=100, it gains $f_{1}$ score of 100%; in its
worst-case scenario, where $\alpha=0.25$ and SNR=50, it respectively gets
98.02% and 99.77% of $f_{1}$ score and accuracy. In the following, the paper
conducts a comparison to figure out how these numbers are meaningful and
whether the proposed model can better mitigate the adverse impacts of
imbalanced and noisy conditions.

###  5.4 Model performance evaluation

In order to achieve meaningful comparisons, some novel FDD frameworks were
employed to perform the diagnosis at different scenarios. CLSTM, df-CNN, sdAE,
GAN-CNN, WELM, CNN-STFT, and CNN-FFT have shown promising performances in the
literature, hence, they were selected for this purpose. Three traditional
machine learning classifiers, SVM, ANN and Random Forest (RF), are also
considered in this experimental comparison. A grid-search was designed on the
hyper-parameters of these models to achieve higher performances. Specifically,
the learning rate, batch size and the architecture of fully connected layers
were optimized for each algorithm, while the number of epochs was set to 50
for all architectures. CLSTM-ELM was also added to the comparison panel to
examine the necessity of Weighted ELM in the architecture of the proposed
framework. This paper used two augmentation techniques: (i) "GAN": with WGAN-
GP, as discussed earlier, (ii) "classic": where the samples are flipped,
mirrored and different white noises are added to the samples. A brief
description of the selected frameworks is provided in Table 2.

Table 2: The comparison panel

Framework |  Augmentation |  Preprossecing |  Description |  Reference
---|---|---|---|---
CLSTM |  classic |  FFT+CWT+Statistical features |  Its architecture comprises two CNN blocks (containing 1D$\sim$-Convolutional layers, Batch Normalization, ReLU and Max Pooling), a LSTM block, a Logarithmic SoftMax, a concatenation which adds statistical features and three fully connected neural networks for the classification. |  [46]
CLSTM-ELM |  classic |  FFT+CWT+Statistical features |  Its CNN and LSTM architecture are the same as in CLSTM; yet, the fully connected layers are substituted for W-ELM with 150 nodes. |  N/A
df-CNN |  classic |  raw signals |  It is proposed to make an abstract 2-dimensional image out of raw signals. Its architecture comprises two CNN blocks (containing 2D$\sim$-Convolutional layers, Batch Normalization, ReLU and Max Pooling), and three fully connected neural networks for the classification. df-CNN works directly on the raw vibration signals. |  [48]
sdAE |  classic |  raw signals |  It is a multilayered architecture composed of four auto-associative neural network layers, which contain one input layer and three AEs. The input of this framework are raw signals. |  [14]
CNN-STFT |  classic |  STFT |  The architecture consists of three CNN blocks (containing one 1D-$\sim$Convolutional layer and a Pooling layer), two fully connected layers, and a SoftMax classification layer. It takes short-term Fourier transform (STFT) form of the signals as its input. |  [31]
GAN-CNN |  WGAN-GP |  STFT |  The paper has slightly modified the architecture proposed by its authors to examine how the same WGAN-GP architecture improves the CNN-STFT model. Therefore, GAN-CNN is a combination of WGAN-GP and CNN-STFT models. |  [33]
W-ELM |  |  FFT+VMD+Statistical features |  It takes a combination of FFT, VMD [49] and some statistical features. |  [20]
CNN-FFT |  classic |  FFT |  It takes FFT as its input while its architecture comprises two CNN blocks (containing 1D-$\sim$Convolutional and Pooling layers) followed by three fully connected layers. |  [16]
GAN-SVM |  WGAN-GP |  Statistical features |  SVM with polynomial kernel and degree of 2 is selected |  N/A
GAN-ANN |  WGAN-GP |  Statistical features |  3 fully connected layers with a grid search to find optimal number of neurons per layer and the activation functions |  N/A
GAN-RF |  WGAN-GP |  Statistical features |  A grid search is designed to find the optimal number of estimators, and criteria (between ‘gini’ and ‘entropy’) parameters |  N/A

To avoid the weight initialization effect and randomness on the results, we
ran each framework for ten independent times, using a k-fold cross validation
technique on each imbalance and noise degree conditions. Figure 11 illustrates
the corresponding normalized confusion matrices and the model performances.
Comparing the different scenarios, it can plainly be concluded that GAN-CLSTM-
ELM has a better ability to extenuate the negative effects of imbalance and
noise conditions compared to the other frameworks. Regarding the highly
imbalanced situation, its $f_{1}$ score has gently dropped by 0.32% in the
first two scenarios (SNR:100, $\alpha:2^{-2}$ and SNR:100, $\alpha:2^{2}$)
while the other frameworks have shown relatively substantial declines in their
$f_{1}$ scores, ranging from 1.14% (GAN-CNN) to roughly 48% (df-CNN). In the
second scenario, while the proposed model correctly identifies all the
minority class samples, CLSTM-ELM and GAN-CNN were able to classify roughly
92% of them. This percentage for CLSTM, sdAE, WELM and CNN-STFT was between 80
and 85. CNN-FFT and df-CNN were among the poorest models to identify the
minority class as df-CNN could not correctly diagnose any of the corresponding
samples. The figure also shows that replacing the fully connected layers with
W-ELM in the CLSTM-ELM model has slightly increased its robustness when
$\alpha$ plummets from 4 to 0.25.

Figure 11: Confusion matrices and $f_{1}$-scores of the comparison panel in
different scenarios(t represents true labeled samples)

In the presence of heavy noises, there are sudden falls in the performances of
all the algorithms. Comparing the first and the third scenarios (SNR=100,
$\alpha=2^{2}$ and SNR=10, $\alpha=2^{2}$), GAN-CLSTM-ELM, CLSTM-ELM and GAN-
CNN had the least decrease in the $f_{1}$ score (roughly 5%); thus, they were
the most robust algorithms in noisy conditions. However, CNN-STFT achieved
comparatively poorer results when $\alpha$ dips below 1. Its combination with
a WGAN-GP mitigated this loss and GAN-CNN achieved a satisfactory result. With
the presence of heavy noises, WELM classification quality drastically plunged
and, despite its comparatively satisfactory performance in the first two
scenarios, the noise made it unable to diagnose the minority class in highly
imbalanced situations. By comparing the confusion matrices of WELM and CLSTM-
ELM, it was demonstrated that CLSTM architecture alongside WELM model improves
its performance against the noise. It is worth noting that both models that
comprise WGAN-GP showed the highest scores in the first, second and the last
scenario and exhibited superiority over their root algorithms, CLSTM-ELM,
CLSTM, WELM and CNN-STFT. This proves that WGAN-GP can effectively enhance the
quality of the classifier not only in imbalanced situations but also in noisy
environments.

##  6 Discussion and Conclusions

In many real applications of fault detection and diagnosis data tend to be
imbalanced and noisy, meaning that the number of samples for some fault
classes is much fewer than the normal data samples and there are errors in
recording the actual measurement by the sensors. These two conditions make
many traditional FDD frameworks perform poorly in real-world industrial
environments.

In this paper a novel framework called GAN-CLSTM-ELM is proposed, which
enhances the performance of rotating machinery FDD systems coping with highly-
imbalanced and noisy datasets. In this framework, WGAN-GP is first applied to
augment the minority class and enhance the training set. A hybrid classifier
is then developed, containing Convolutional LSTM and Weighted ELM, which
learns more efficiently from vibration signals. The framework also benefits
from both wavelet and Fourier transform techniques in its feature engineering
step, revealing more hidden information of the fault signatures to make the
classifier perform more accurately. The effectiveness of the proposed
framework is verified by using four dataset settings with different imbalance
severities and SNRs. After conducting the comparisons with state-of-the-art
FDD algorithms, it is demonstrated that the GAN-CLSTM-ELM framework can reduce
the misclassification rate and outperform the other methods, more
significantly when the imbalance degree is higher. The efficiency of the WGAN-
GP is also proved by comparing the results of the proposed model and CLSTM-ELM
as well as those of GAN-CNN and CNN-SFTF models. The experimental results make
it discernible that using a generative algorithm helps to alleviate the
adverse impacts of low SNRs. Therefore, it stresses the necessity of employing
such hybrid frameworks for practitioners working on noisy and industrial
applications. The paper also justifies the implementation of W-ELM in the
architecture of CLSTM, since the adjusted model shows sturdy classification
when $\alpha$ decreases either in noisy or noiseless scenarios. A sensitivity
analysis is designed with 25 dataset settings built on a range of $\alpha$ and
SNR values, to obtain insights of how these two factors impact on the model’s
classification ability.

Extracting the FFT and CWT spectra needs some knowledge of signal processing
and is still more convenient than extracting other hand-crafted features
proposed in the literature. Another advantage of the proposed framework is
that it gains comparatively high performances under noisy conditions while it
requires no complex denoising pre-processing being handled by employees with
expert knowledge of signal processing. These characteristics make GAN-CLSTM-
ELM an attractive option for industrial practitioners who are in need of a
relatively easy-to-use software without undergoing any complicated pre-
processing task.

Future work will include more experiments on the behavior of different
generative algorithms and the development of a more powerful architecture to
create high-quality signals with fewer samples. We will also attempt to
explore the feasibility of implementing and testing the proposed framework on
other applications.

## References

  * Glowacz et al. [2018] A. Glowacz, W. Glowacz, Z. Glowacz, J. Kozik,  Early fault diagnosis of bearing and stator faults of the single-phase induction motor using acoustic signals,  Measurement: Journal of the International Measurement Confederation 113 (2018) 1–9. doi:10.1016/j.measurement.2017.08.036.
  * Shojaeinasab et al. [2022] A. Shojaeinasab, T. Charter, M. Jalayer, M. Khadivi, O. Ogunfowora, N. Raiyani, M. Yaghoubi, H. Najjaran,  Intelligent manufacturing execution systems: A systematic review,  Journal of Manufacturing Systems 62 (2022) 503–522. URL: https://doi.org/10.1016/j.jmsy.2022.01.004. doi:10.1016/j.jmsy.2022.01.004.
  * Liu et al. [2018] R. Liu, B. Yang, E. Zio, X. Chen,  Artificial intelligence for fault diagnosis of rotating machinery: A review,  Mechanical Systems and Signal Processing 108 (2018) 33–47. URL: https://doi.org/10.1016/j.ymssp.2018.02.016. doi:10.1016/j.ymssp.2018.02.016.
  * Zhao et al. [2016] G. Zhao, G. Zhang, Q. Ge, X. Liu,  Research Advances in Fault Diagnosis and Prognostic based on Deep Learning,  Prognostics and System Health Management Conference (2016) 1–6. doi:10.1109/PHM.2016.7819786.
  * Balderston [1969] H. Balderston, Incipient failure detection: Incipient failure detection in ball bearings, Technical Report, BOEING CO SEATTLE WA AEROSPACE SYSTEMS DIV., 1969.
  * Weichbrodt and Smith [1970] B. Weichbrodt, K. A. Smith,  Signature Analaysis. Non-intrusive techniques for incipient failure identification,  in: NBS Spec Publ 336, Proc 5th Space Simulation Symp Conf, Elsevier B.V., Gaithersburg, MD, USA, 1970.
  * Li and Ma [1997] C. J. Li, J. Ma,  Wavelet decomposition of vibrations for detection of bearing-localized defects 30 (1997) 143–149. doi:10.1016/S0963-8695(96)00052-7.
  * Raghav and Sharma [2020] M. S. Raghav, R. B. Sharma,  A Review on Fault Diagnosis and Condition Monitoring of Gearboxes by Using AE Technique,  Archives of Computational Methods in Engineering (2020). doi:10.1007/s11831-020-09480-8.
  * Mao et al. [2017] W. Mao, L. He, Y. Yan, J. Wang,  Online sequential prediction of bearings imbalanced fault diagnosis by extreme learning machine,  Mechanical Systems and Signal Processing 83 (2017) 450–473. URL: http://dx.doi.org/10.1016/j.ymssp.2016.06.024. doi:10.1016/j.ymssp.2016.06.024.
  * Lee [2017] H. Lee,  Framework and development of fault detection classification using IoT device and cloud environment,  Journal of Manufacturing Systems 43 (2017) 257–270. URL: http://dx.doi.org/10.1016/j.jmsy.2017.02.007. doi:10.1016/j.jmsy.2017.02.007.
  * Zhang et al. [2016] T. Zhang, H. Ye, H. Zhang, M. Li,  PCA-LMNN-Based Fault Diagnosis Method for Ironmaking Processes with Insufficient Faulty Data,  ISIJ International 56 (2016) 1779–1788. URL: https://www.jstage.jst.go.jp/article/isijinternational/56/10/56{_}ISIJINT-2016-101/{_}article. doi:10.2355/isijinternational.ISIJINT-2016-101.
  * Liu et al. [2018] G. Liu, H. Bao, B. Han,  A Stacked Autoencoder-Based Deep Neural Network for Achieving Gearbox Fault Diagnosis,  Mathematical Problems in Engineering 2018 (2018) 1–10. URL: https://www.hindawi.com/journals/mpe/2018/5105709/. doi:10.1155/2018/5105709.
  * Haidong et al. [2017] S. Haidong, J. Hongkai, Z. Huiwei, W. Fuan,  A novel deep autoencoder feature learning method for rotating machinery fault diagnosis,  Mechanical Systems and Signal Processing 95 (2017) 187–204. URL: http://dx.doi.org/10.1016/j.ymssp.2017.03.034. doi:10.1016/j.ymssp.2017.03.034.
  * Lu et al. [2017] C. Lu, Z. Y. Wang, W. L. Qin, J. Ma,  Fault diagnosis of rotary machinery components using a stacked denoising autoencoder-based health state identification,  Signal Processing 130 (2017) 377–388. URL: http://dx.doi.org/10.1016/j.sigpro.2016.07.028. doi:10.1016/j.sigpro.2016.07.028.
  * Ke Li and Quanxin Wang [2015] Ke Li, Quanxin Wang,  Study on signal recognition and diagnosis for spacecraft based on deep learning method,  in: 2015 Prognostics and System Health Management Conference (PHM), IEEE, 2015, pp. 1–5. URL: http://ieeexplore.ieee.org/document/7380040/. doi:10.1109/PHM.2015.7380040.
  * Jia et al. [2018] F. Jia, Y. Lei, L. Guo, J. Lin, S. Xing,  A neural network constructed by deep learning technique and its application to intelligent fault diagnosis of machines,  Neurocomputing 272 (2018) 619–628. URL: https://doi.org/10.1016/j.neucom.2017.07.032. doi:10.1016/j.neucom.2017.07.032.
  * Lei et al. [2016] Y. Lei, F. Jia, J. Lin, S. Xing, S. X. Ding,  An Intelligent Fault Diagnosis Method Using Unsupervised Feature Learning Towards Mechanical Big Data,  IEEE Transactions on Industrial Electronics 63 (2016) 3137–3147. URL: http://ieeexplore.ieee.org/document/7386639/. doi:10.1109/TIE.2016.2519325.
  * Li et al. [2014] K. Li, X. Kong, Z. Lu, L. Wenyin, J. Yin,  Boosting weighted ELM for imbalanced learning,  Neurocomputing 128 (2014) 15–21. URL: http://dx.doi.org/10.1016/j.neucom.2013.05.051. doi:10.1016/j.neucom.2013.05.051.
  * Vong et al. [2014] C. M. Vong, W. F. Ip, P. K. Wong, C. C. Chiu,  Predicting minority class for suspended particulate matters level by extreme learning machine,  Neurocomputing 128 (2014) 136–144. URL: http://dx.doi.org/10.1016/j.neucom.2012.11.056. doi:10.1016/j.neucom.2012.11.056.
  * Zong et al. [2013] W. Zong, G. B. Huang, Y. Chen,  Weighted extreme learning machine for imbalance learning,  Neurocomputing 101 (2013) 229–242. URL: http://dx.doi.org/10.1016/j.neucom.2012.08.010. doi:10.1016/j.neucom.2012.08.010.
  * Mirza et al. [2013] B. Mirza, Z. Lin, K. A. Toh,  Weighted online sequential extreme learning machine for class imbalance learning,  Neural Processing Letters 38 (2013) 465–486. doi:10.1007/s11063-013-9286-9.
  * Mirza et al. [2015] B. Mirza, Z. Lin, J. Cao, X. Lai,  Voting based weighted online sequential extreme learning machine for imbalance multi-class classification,  Proceedings - IEEE International Symposium on Circuits and Systems 2015-July (2015) 565–568. doi:10.1109/ISCAS.2015.7168696.
  * Hao and Liu [2020] W. Hao, F. Liu,  Imbalanced Data Fault Diagnosis Based on an Evolutionary Online Sequential Extreme Learning Machine,  Symmetry 12 (2020) 1204. URL: https://www.mdpi.com/2073-8994/12/8/1204. doi:10.3390/sym12081204.
  * Chen et al. [2017] Y.-J. Chen, B.-c. Wang, J.-Z. Wu, C.-F. Chien,  Big Data Analytic for Multivariate Fault Detection and Classification in Semiconductor Manufacturing,  2017\.
  * Cheng et al. [2019] J. Cheng, J. Chen, Y. nan Guo, S. Cheng, L. Yang, P. Zhang,  Adaptive CCR-ELM with variable-length brain storm optimization algorithm for class-imbalance learning,  Natural Computing 9 (2019). URL: https://doi.org/10.1007/s11047-019-09735-9. doi:10.1007/s11047-019-09735-9.
  * Xu et al. [2020] Q. Xu, S. Lu, W. Jia, C. Jiang,  Imbalanced fault diagnosis of rotating machinery via multi-domain feature extraction and cost-sensitive learning,  Journal of Intelligent Manufacturing 31 (2020) 1467–1481. URL: https://doi.org/10.1007/s10845-019-01522-8. doi:10.1007/s10845-019-01522-8.
  * Zhao et al. [2020] X. Zhao, M. Jia, Z. Liu,  Fault Diagnosis Framework of Rolling Bearing Using Adaptive Sparse Contrative Auto-Encoder with Optimized Unsupervised Extreme Learning Machine,  IEEE Access 8 (2020) 99154–99170. doi:10.1109/ACCESS.2019.2963193.
  * Zhao et al. [2019] Y. P. Zhao, G. Huang, Q. K. Hu, J. F. Tan, J. J. Wang, Z. Yang,  Soft extreme learning machine for fault detection of aircraft engine,  Aerospace Science and Technology 91 (2019) 70–81. URL: https://doi.org/10.1016/j.ast.2019.05.021. doi:10.1016/j.ast.2019.05.021.
  * Han et al. [2020] B. Han, S. Jia, G. Liu, J. Wang,  Imbalanced Fault Classification of Bearing via Wasserstein Generative Adversarial Networks with Gradient Penalty,  Shock and Vibration 2020 (2020) 1–14. URL: https://www.hindawi.com/journals/sv/2020/8836477/. doi:10.1155/2020/8836477.
  * Ding et al. [2019] Y. Ding, L. Ma, J. Ma, C. Wang, C. Lu,  A generative adversarial network-based intelligent fault diagnosis method for rotating machinery under small sample size conditions,  IEEE Access 7 (2019) 149736–149749. doi:10.1109/ACCESS.2019.2947194.
  * Zhang et al. [2020] Y. Zhang, K. Xing, R. Bai, D. Sun, Z. Meng,  An enhanced convolutional neural network for bearing fault diagnosis based on time–frequency image,  Measurement: Journal of the International Measurement Confederation 157 (2020) 107667\. URL: https://doi.org/10.1016/j.measurement.2020.107667. doi:10.1016/j.measurement.2020.107667.
  * He et al. [2020] W. He, Y. He, B. Li,  Generative Adversarial Networks with Comprehensive Wavelet Feature for Fault Diagnosis of Analog Circuits,  IEEE Transactions on Instrumentation and Measurement 69 (2020) 6640–6650. doi:10.1109/TIM.2020.2969008.
  * Liang et al. [2020] P. Liang, C. Deng, J. Wu, Z. Yang,  Intelligent fault diagnosis of rotating machinery via wavelet transform, generative adversarial nets and convolutional neural network,  Measurement: Journal of the International Measurement Confederation 159 (2020) 107768\. URL: https://doi.org/10.1016/j.measurement.2020.107768. doi:10.1016/j.measurement.2020.107768.
  * Goodfellow et al. [2014] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, Y. Bengio,  Generative adversarial nets,  Advances in neural information processing systems 27 (2014).
  * Ren et al. [2019] J. Ren, Y. Liu, J. Liu,  EWGAN: Entropy-Based Wasserstein GAN for Imbalanced Learning,  Proceedings of the AAAI Conference on Artificial Intelligence 33 (2019) 10011–10012. doi:10.1609/aaai.v33i01.330110011.
  * Arjovsky et al. [2017] M. Arjovsky, S. Chintala, L. Bottou,  Wasserstein GaN,  arXiv (2017). arXiv:1701.07875.
  * Zhou et al. [2019] F. Zhou, X. Lin, C. Liu, Y. Zhao, P. Xu, L. Ren, T. Xue, L. Ren,  A survey of visualization for smart manufacturing,  Journal of Visualization 22 (2019) 419–435. URL: https://doi.org/10.1007/s12650-018-0530-2. doi:10.1007/s12650-018-0530-2.
  * Jalayer et al. [2021] M. Jalayer, R. Jalayer, A. Kaboli, C. Orsenigo, C. Vercellis,  Automatic visual inspection of rare defects: A framework based on gp-wgan and enhanced faster r-cnn,  in: 2021 IEEE International Conference on Industry 4.0, Artificial Intelligence, and Communications Technology (IAICT), IEEE, 2021, pp. 221–227. URL: https://doi.org/10.1109/IAICT52856.2021.9532584. doi:10.1109/IAICT52856.2021.9532584.
  * Gers et al. [1999] F. A. Gers, J. Schmidhuber, F. Cummins,  Learning to forget: Continual prediction with LSTM,  IEE Conference Publication 2 (1999) 850–855. doi:10.1049/cp:19991218.
  * Gers and Schraudolph [2002] F. A. Gers, N. N. Schraudolph,  Learning Precise Timing with LSTM Recurrent Networks,  Journal of Machine Learning Research 3 (2002) 115–143.
  * Healy et al. [1972] M. J. R. Healy, C. R. Rao, S. K. Mitra,  Generalized Inverse of Matrices and its Applications.,  Journal of the Royal Statistical Society. Series A (General) 135 (1972) 439. URL: https://www.jstor.org/stable/10.2307/2344631?origin=crossref. doi:10.2307/2344631.
  * Hoerl and Kennard [1970] A. E. Hoerl, R. W. Kennard,  Ridge Regression: Biased Estimation for Nonorthogonal Problems,  Technometrics 12 (1970) 55–67. URL: http://www.tandfonline.com/doi/abs/10.1080/00401706.1970.10488634. doi:10.1080/00401706.1970.10488634.
  * Fletcher [2000] R. Fletcher, Practical Methods of Optimization, Second Edition, John Wiley & Sons, Ltd, Chichester, West Sussex England, 2000\. URL: http://doi.wiley.com/10.1002/9781118723203. doi:10.1002/9781118723203.
  * Gulrajani et al. [2017] I. Gulrajani, F. Ahmed, M. Arjovsky, V. Dumoulin, A. Courville,  Improved training of wasserstein GANs,  Advances in Neural Information Processing Systems 2017-Decem (2017) 5768–5778. arXiv:1704.00028.
  * Mallat [2008] S. Mallat, A Wavelet Tour of Signal Processing, 3rd editio ed., Academic Press, 2008. URL: https://www.elsevier.com/books/a-wavelet-tour-of-signal-processing/mallat/978-0-12-374370-1.
  * Jalayer et al. [2021] M. Jalayer, C. Orsenigo, C. Vercellis,  Fault detection and diagnosis for rotating machinery: A model based on convolutional LSTM, Fast Fourier and continuous wavelet transforms,  Computers in Industry 125 (2021) 103378. URL: https://doi.org/10.1016/j.compind.2020.103378. doi:10.1016/j.compind.2020.103378.
  * Karim et al. [2019] F. Karim, S. Majumdar, H. Darabi,  Insights into lstm fully convolutional networks for time series classification,  IEEE Access 7 (2019) 67718–67725. doi:10.1109/ACCESS.2019.2916828. arXiv:1902.10756.
  * Mao et al. [2021] W. Mao, W. Feng, Y. Liu, D. Zhang, X. Liang,  A new deep auto-encoder method with fusing discriminant information for bearing fault diagnosis,  Mechanical Systems and Signal Processing 150 (2021) 107233. URL: https://doi.org/10.1016/j.ymssp.2020.107233. doi:10.1016/j.ymssp.2020.107233.
  * Yan and Jia [2018] X. Yan, M. Jia,  A novel optimized SVM classification algorithm with multi-domain feature and its application to fault diagnosis of rolling bearing,  Neurocomputing 313 (2018) 47–64. URL: https://doi.org/10.1016/j.neucom.2018.05.002. doi:10.1016/j.neucom.2018.05.002.


$\frac{dx}{dt}$

DIFFERENTIAL EQUATIONS

AND

CONTROL PROCESSES

N. 4, 2020

Electronic Journal,
reg. N ${\Phi}$C77-39410 at 15.04.2010

ISSN 1817-2172

http://diffjournal.spbu.ru/

e-mail<EMAIL_ADDRESS>
Numerical methods

# A Numerical-Analytical Method for Constructing Periodic Solutions of the
Lorenz System

Alexander N. Pchelintsev
Tambov State Technical University,
ul. Sovetskaya 106, Tambov, 392000, Russia
e-mail<EMAIL_ADDRESS>
Abstract. This article describes a method for constructing approximations to
periodic solutions of dynamic Lorenz system with classical values of the
system parameters. The author obtained a system of nonlinear algebraic
equations in general form concerning of the cyclic frequency, constant terms
and amplitudes of harmonics that make up harmonic approximations to the
desired solutions. The initial approximation for the Newton method is
selected, which converges to a solution describing a periodic solution
different from the equilibrium position. The results of a computational
experiment are presented. The results are verified using high-precision
calculations.

Keywords: Attractor, Lorenz Attractor, Trigonometric Polynomial, Newton’s
Method.

##  1 Introduction

Let us consider the nonlinear system of differential equations introduced by
E. Lorenz in [1]

$\left\\{\begin{array}[]{l}\dot{x}_{1}=\sigma(x_{2}-x_{1}),\vskip 6.0pt plus
2.0pt minus 2.0pt\\\ \dot{x}_{2}=rx_{1}-x_{2}-x_{1}x_{3},\vskip 6.0pt plus
2.0pt minus 2.0pt\\\ \dot{x}_{3}=x_{1}x_{2}-bx_{3},\end{array}\right.$ (1)

where $\sigma=10$, $r=28$ and $b=8/3$ are the classical values of the system
parameters.

Let us denote by
$X(t)=\left[x_{1}(t)\>\>x_{2}(t)\>\>x_{3}(t)\right]^{\scriptsize\mbox{T}}$. It
is proved in the article [1] that there exists a number $C>0$ such that for
any solution $X(t)$ of the system (1), starting at time moment, $|X(t)|<C$,
and the divergence of the vector velocity field of the system (1) is negative
everywhere in $\mathbb{R}^{3}$ for classical values of the system parameters.
Then [1] there exists a limit set, called the Lorenz attractor, to which all
trajectories of the dynamical system are attracted when time tends to
infinity. Thus the attractor determines the behavior of the solutions of a
dynamical system over large segments of time.

W. Tucker in his work [2] proved that the attractor is hyperbolic in the
system (1), that is, the attractor consists of cycles everywhere dense on it
along which the near trajectories diverge exponentially. This creates their
chaotic behavior.

As know [3, 4], the symbolic dynamics is used to track cycles in the Lorenz
system. The region in the phase space containing the attractor is divided into
a finite number of subdomains. Denoting each partition element by a symbol,
the trajectories on the attractor passing through the corresponding regions
are encoded by sequences of such symbols. If the sequence has regularity
(repeatability of groups of characters), then we have a cycle. However, the
return of trajectories in a neighborhood of its part does not mean its
closure. A critique of the results of such computational experiments can be
found, for example, in [5].

In 2004, D. Viswanath published the paper [6], in which he presented the
initial conditions and periods for three cycles in the Lorenz attractor with a
high accuracy. The calculation algorithm is based on the Lindstedt-Poincaré
(LP) method, which (unlike numerical integration methods) is not affected by
the stability of the cycle to which approximations are constructed.

An analysis of the Viswanath’s articles [6, 7] showed that the author gives a
general description of the algorithm without reference to the computer
implementation (in MATLAB as indicated in his works). Moreover, it is not
clear how the obtained inhomogeneous linear system of differential equations
with periodic coefficients is symbolically solved by the LP-method. For
example, this can be done for the Van der Pol equation without any special
problems.

In the article [6] Viswanath showed data that can be verified by solving the
Cauchy problem with high-precision numerical methods (for example, [8]), but
the details of the algorithm are not disclosed.

Therefore, it is important here to obtain the values of the initial conditions
and the period with a given accuracy, having described in detail the
implementation of the cycles search algorithm in the system (1).

The goal of this article is to develop a numerical-analytical method for
constructing approximations to periodic solutions of the Lorenz system, which
is simpler to implement than the LP-method. In this case, a system of
nonlinear algebraic equations concerning of the cyclic frequency, constant
terms, and amplitudes of harmonics making up the desired solution will be
obtained in general form.

##  2 A Numerical-Analytical Method

Attempts to construct approximate periodic solutions in the system (1) with
were made before Viswanath (for example, [9]) by the method of harmonic
balance, but with low accuracy in representing real numbers, while in the
article [9] initial conditions and periods of found cycles are not indicated
(only drawings with cycles are given). Now this method is actively developing
in the works of [10, 11, 12] A. Luo to find periodic solutions of nonlinear
systems of differential equations.

Next, we describe a numerical-analytical method for constructing
approximations to periodic solutions of the system (1). We make for this an
approximation of the phase coordinates on the period $T$ by trigonometric
polynomials in general form with an unknown cyclic frequency $\omega$ (since
we do not know the value of $T$; in the general case, it can be an irrational
number):

$\begin{array}[]{l}\displaystyle
x_{1}(t)\approx\tilde{x}_{1}(t)=x_{1,0}+\sum_{i=1}^{h}\left(c_{1,i}\cos(i\omega
t)+s_{1,i}\sin(i\omega t)\right),\vskip 6.0pt plus 2.0pt minus 2.0pt\\\
\displaystyle
x_{2}(t)\approx\tilde{x}_{2}(t)=x_{2,0}+\sum_{i=1}^{h}\left(c_{2,i}\cos(i\omega
t)+s_{2,i}\sin(i\omega t)\right),\vskip 6.0pt plus 2.0pt minus 2.0pt\\\
\displaystyle
x_{3}(t)\approx\tilde{x}_{3}(t)=x_{3,0}+\sum_{i=1}^{h}\left(c_{3,i}\cos(i\omega
t)+s_{3,i}\sin(i\omega t)\right),\end{array}$

where $h$ is given number of harmonics. If $i>h$, then we assume

$c_{1,i}=s_{1,i}=c_{2,i}=s_{2,i}=c_{3,i}=s_{3,i}=0.$ (2)

By the right-hand side of the system (1), we compose the residuals

$\begin{array}[]{l}\delta_{1}(t)=\tilde{x}^{\prime}_{1}(t)-\sigma[\tilde{x}_{2}(t)-\tilde{x}_{1}(t)],\\\
\delta_{2}(t)=\tilde{x}^{\prime}_{2}(t)-[r\tilde{x}_{1}(t)-\tilde{x}_{2}(t)-\tilde{x}_{1}(t)\tilde{x}_{3}(t)],\\\
\delta_{3}(t)=\tilde{x}^{\prime}_{3}(t)-[\tilde{x}_{1}(t)\tilde{x}_{2}(t)-b\tilde{x}_{3}(t)],\end{array}$

where the prime denotes the time derivative of the function. If we make
calculations in an analytical form, then for each residual you need the
following:

  1. 1.

Differentiate by time the corresponding trigonometric polynomial;

  2. 2.

Where there are products of phase coordinates, multiply the corresponding
trigonometric polynomials, converting the products of trigonometric functions
into sums;

  3. 3.

Give similar terms for each function $\cos()$ and $\sin()$ with the
corresponding argument;

  4. 4.

By virtue of the equalities (2), to cut off the higher-order harmonics from
the resulting residual;

  5. 5.

Set the resulting residual to zero, i.e., coefficients at its harmonics.

If we put together the found algebraic equations for each residual, we obtain
a still unclosed system of nonlinear equations concerning of unknown
amplitudes $c_{1,i}$, $s_{1,i}$, $c_{2,i}$, $s_{2,i}$, $c_{3,i}$ and $s_{3,i}$
($i=\overline{1,h}$), constant terms $x_{1,0}$, $x_{2,0}$ and $x_{3,0}$ and
the cyclic frequency $\omega$. The number of unknown variables in the system
is $3(1+2h)+1=6h+4$, but the equations are one less.

An additional equation can be taken from the following considerations. It is
known (see [4, 6]) that the desired cycles intersect the plane passing through
the equilibrium positions of the system (1)

$O_{1}\left(-\sqrt{b(r-1)},\,-\sqrt{b(r-1)},\,r-1\right),\>\>O_{2}\left(\sqrt{b(r-1)},\,\sqrt{b(r-1)},\,r-1\right)$
(3)

and parallel to the plane $x_{1}Ox_{2}$ (a Poincare section). Then the third
coordinate in the initial condition for the desired cycles is equal to $r-1$,
whence $\tilde{x}_{3}(0)=r-1$.

Therefore the additional equation of the system has the form:

$x_{3,0}+\sum_{i=1}^{h}c_{3,i}-27=0.$

The author did not find in literature of other additional information on the
periodic solutions in the Lorenz system. Note that for the three cycles found
by Viswanath, in the initial condition for the third coordinate, the number 27
was taken.

Next, we give an example of a system of equations for $h=2$:

$\left\\{\begin{aligned} \omega s_{1,1}-10c_{2,1}+10c_{1,1}&=0,\\\
-10s_{2,1}+10s_{1,1}-c_{1,1}\omega&=0,\\\ 2\omega
s_{1,2}-10c_{2,2}+10c_{1,2}&=0,\\\ -10s_{2,2}+10s_{1,2}-2c_{1,2}\omega&=0,\\\
10x_{1,0}-10x_{2,0}&=0,\\\
c_{1,1}x_{3,0}+c_{3,1}x_{1,0}+\dfrac{s_{1,1}s_{3,2}}{2}+\dfrac{s_{1,2}s_{3,1}}{2}+\omega
s_{2,1}+\dfrac{c_{1,1}c_{3,2}}{2}+\dfrac{c_{1,2}c_{3,1}}{2}+c_{2,1}-28c_{1,1}&=0,\\\
s_{1,1}x_{3,0}+s_{3,1}x_{1,0}+\dfrac{c_{1,1}s_{3,2}}{2}-\dfrac{c_{1,2}s_{3,1}}{2}+s_{2,1}+\dfrac{c_{3,1}s_{1,2}}{2}-\dfrac{c_{3,2}s_{1,1}}{2}-28s_{1,1}-c_{2,1}\omega&=0,\\\
c_{1,2}x_{3,0}+c_{3,2}x_{1,0}-\dfrac{s_{1,1}s_{3,1}}{2}+2\omega
s_{2,2}+\dfrac{c_{1,1}c_{3,1}}{2}+c_{2,2}-28c_{1,2}&=0,\\\
s_{1,2}x_{3,0}+s_{3,2}x_{1,0}+\dfrac{c_{1,1}s_{3,1}}{2}+s_{2,2}-28s_{1,2}+\dfrac{c_{3,1}s_{1,1}}{2}-2c_{2,2}\omega&=0,\\\
x_{1,0}x_{3,0}+x_{2,0}-28x_{1,0}+\dfrac{s_{1,2}s_{3,2}}{2}+\dfrac{s_{1,1}s_{3,1}}{2}+\dfrac{c_{1,2}c_{3,2}}{2}+\dfrac{c_{1,1}c_{3,1}}{2}&=0,\\\
-c_{1,1}x_{2,0}-c_{2,1}x_{1,0}+\omega
s_{3,1}-\dfrac{s_{1,1}s_{2,2}}{2}-\dfrac{s_{1,2}s_{2,1}}{2}+\dfrac{8c_{3,1}}{3}-\dfrac{c_{1,1}c_{2,2}}{2}-\dfrac{c_{1,2}c_{2,1}}{2}&=0,\\\
-s_{1,1}x_{2,0}-s_{2,1}x_{1,0}+\dfrac{8s_{3,1}}{3}-\dfrac{c_{1,1}s_{2,2}}{2}+\dfrac{c_{1,2}s_{2,1}}{2}-\dfrac{c_{2,1}s_{1,2}}{2}+\dfrac{c_{2,2}s_{1,1}}{2}-c_{3,1}\omega&=0,\\\
-c_{1,2}x_{2,0}-c_{2,2}x_{1,0}+2\omega
s_{3,2}+\dfrac{s_{1,1}s_{2,1}}{2}+\dfrac{8c_{3,2}}{3}-\dfrac{c_{1,1}c_{2,1}}{2}&=0,\\\
-s_{1,2}x_{2,0}-s_{2,2}x_{1,0}+\dfrac{8s_{3,2}}{3}-\dfrac{c_{1,1}s_{2,1}}{2}-\dfrac{c_{2,1}s_{1,1}}{2}-2c_{3,2}\omega&=0,\\\
\dfrac{8x_{3,0}}{3}-x_{1,0}x_{2,0}-\dfrac{s_{1,2}s_{2,2}}{2}-\dfrac{s_{1,1}s_{2,1}}{2}-\dfrac{c_{1,2}c_{2,2}}{2}-\dfrac{c_{1,1}c_{2,1}}{2}&=0,\\\
x_{3,0}+c_{3,1}+c_{3,2}-27&=0.\end{aligned}\right.$

Note that for any $h$ a similar system has solutions

$\begin{array}[]{c}\displaystyle
x_{1,0}=x_{2,0}=\pm\sqrt{b(r-1)},\>x_{3,0}=r-1,\>c_{k,i}=0,\>s_{k,i}=0,\\\
\omega\>\mbox{is any
number},\>\,k=\overline{1,3},\>i=\overline{1,h},\end{array}$

corresponding to the equilibrium positions (3).

Therefore the resulting nonlinear system of algebraic equations has a non-
unique solution. To find its approximate solutions, we will use the Newton
numerical method, whose a convergence to the desired solution (i.e.,
describing a periodic solution of the system (1) different from its the
equilibrium positions) depends on the choice of the initial approximation.

##  3 The Symbolic Computations to Obtain the System of Algebraic Equations

Thus, to obtain an approximation to the periodic solution, we must obtain a
nonlinear system concerning of unknown decomposition coefficients and
frequencies. As shown in the previous section, even for two harmonics, the
system has a bulky form. Therefore, we consider the algorithm for performing
symbolic calculations to obtain it.

When developing software [13], the Maxima math package (a computer algebra
system) was chosen. The program for obtaining the amplitudes and constant
terms of the residuals for $h=2$ is presented below.

    /* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
    /* [wxMaxima: input   start ] */
    display2d:false$
    x1:x10+c1c1*cos(1*omega*t)+s1c1*sin(1*omega*t)+
    c1c2*cos(2*omega*t)+s1c2*sin(2*omega*t)$
    x2:x20+c2c1*cos(1*omega*t)+s2c1*sin(1*omega*t)+
    c2c2*cos(2*omega*t)+s2c2*sin(2*omega*t)$
    x3:x30+c3c1*cos(1*omega*t)+s3c1*sin(1*omega*t)+
    c3c2*cos(2*omega*t)+s3c2*sin(2*omega*t)$
    assume(omega > 0)$
    delta1:trigreduce(diff(x1,t)-(10*(x2-x1)),t)$
    delta2:trigreduce(diff(x2,t)-(28*x1-x2-x1*x3),t)$
    delta3:trigreduce(diff(x3,t)-(x1*x2-8/3*x3),t)$
    expand(diff(delta1,cos(1*omega*t)));
    expand(diff(delta1,sin(1*omega*t)));
    expand(diff(delta1,cos(2*omega*t)));
    expand(diff(delta1,sin(2*omega*t)));
    expand(integrate(delta1,t,0,2*%pi/omega)*omega/(2*%pi));
    expand(diff(delta2,cos(1*omega*t)));
    expand(diff(delta2,sin(1*omega*t)));
    expand(diff(delta2,cos(2*omega*t)));
    expand(diff(delta2,sin(2*omega*t)));
    expand(integrate(delta2,t,0,2*%pi/omega)*omega/(2*%pi));
    expand(diff(delta3,cos(1*omega*t)));
    expand(diff(delta3,sin(1*omega*t)));
    expand(diff(delta3,cos(2*omega*t)));
    expand(diff(delta3,sin(2*omega*t)));
    expand(integrate(delta3,t,0,2*%pi/omega)*omega/(2*%pi));
    /* [wxMaxima: input   end   ] */

The expression `display2d:false$` turns off multi-line drawing of fractions,
degrees, etc. The sign `$` allows to calculate the result of an expression,
but not display it (instead of `;`). The function `trigreduce(expression,t)`
collapses all products of trigonometric functions concerning of the variable
$t$ in a combination of sums. Differentiation of residuals according to
harmonic functions is necessary to obtain the corresponding amplitudes. The
function `expand(expression)` expands brackets (performs multiplication,
exponentiation, leads similar terms).

To find the constant terms of the residuals, their integration over the period
is applied, i.e. the constant term of the $k$-residual is

$\dfrac{\displaystyle\omega\int_{0}^{\frac{2\pi}{\omega}}\delta_{k}(t)dt}{2\pi}.$

So that during symbolic integration the package does not ask a question about
the sign of the frequency, a command is given `assume(omega > 0)$`.

A file with package commands is generated similarly for any number of $h$
harmonics by a computer program written in C++ [13]. After executing this
program, the package will output symbolic expressions to the console for the
left side of the system of algebraic equations, which will be solved in it by
the Newton method.

Note that the most time-consuming operation here is symbolic integration. For
example, for 120 harmonics, the system formation time is more than 2 days. We
can here parallelize the computational process on three computers, but this
will not have a significant effect. Therefore, a system of algebraic equations
must be formed immediately. Next, we get a general form of this system. Note
that when solving the system of nonlinear equations by the Newton method, the
Jacobi matrix for the left side of the system does not invert. The Maxima
package uses LU decomposition to solve a system of linear equations at each
iteration of the method.

##  4 General Form of the System of Algebraic Equations

Since the right-hand side of the (1) system contains nonlinearities in the
form of products of phase coordinates, let us obtain relations expressing the
coefficients of trigonometric polynomials obtained by multiplying the
approximations $\tilde{x}_{1}(t)\tilde{x}_{3}(t)$ and
$\tilde{x}_{1}(t)\tilde{x}_{2}(t)$.

We consider two functions $f(t)$ and $F(t)$ represented by Fourier series

$\begin{array}[]{c}\displaystyle
f(t)=a_{0}+\sum_{i=1}^{\infty}\left(a_{i}\cos(i\omega t)+b_{i}\sin(i\omega
t)\right),\vskip 6.0pt plus 2.0pt minus 2.0pt\\\ \displaystyle
F(t)=A_{0}+\sum_{i=1}^{\infty}\left(A_{i}\cos(i\omega t)+B_{i}\sin(i\omega
t)\right).\end{array}$

Let

$f(t)F(t)=\alpha_{0}+\sum_{i=1}^{\infty}\left(\alpha_{i}\cos(i\omega
t)+\beta_{i}\sin(i\omega t)\right).$

Following the book [14, pp. 123-125], we have the relations:

$\alpha_{0}=a_{0}A_{0}+\dfrac{1}{2}\sum_{m=1}^{\infty}\left(a_{m}A_{m}+b_{m}B_{m}\right),$
$\alpha_{i}=a_{0}A_{i}+\dfrac{1}{2}\sum_{m=1}^{\infty}\left(a_{m}(A_{m+i}+A_{m-i})+b_{m}(B_{m+i}+B_{m-i})\right),$
(4)
$\beta_{i}=a_{0}B_{i}+\dfrac{1}{2}\sum_{m=1}^{\infty}\left(a_{m}(B_{m+i}-B_{m-i})-b_{m}(A_{m+i}-A_{m-i})\right).$
(5)

We assume that for $i>h$

$a_{i}=b_{i}=A_{i}=B_{i}=0.$

Since for our problem we find for an approximation up to and including the
$h$-harmonic, we zero all the amplitudes in the product for $i>h$, i.e.

$\alpha_{i}=\beta_{i}=0.$

Thus, we pass from the product of series to the product of trigonometric
polynomials. Also in the relations (4) and (5) we will assume [14, p. 124]
that

$A_{m-i}=A_{i-m},\>\>B_{m-i}=-B_{i-m},\>\>B_{0}=0.$

Then we get

$\alpha_{0}=a_{0}A_{0}+\dfrac{1}{2}\sum_{m=1}^{h}\left(a_{m}A_{m}+b_{m}B_{m}\right),$
$\displaystyle\alpha_{i}$
$\displaystyle=a_{0}A_{i}+\dfrac{1}{2}\sum_{m=1}^{\infty}a_{m}A_{m+i}+\dfrac{1}{2}\sum_{m=1}^{\infty}a_{m}A_{m-i}+\dfrac{1}{2}\sum_{m=1}^{\infty}b_{m}B_{m+i}+\dfrac{1}{2}\sum_{m=1}^{\infty}b_{m}B_{m-i}=$
$\displaystyle=a_{0}A_{i}+\dfrac{1}{2}\sum_{m=1}^{h-i}a_{m}A_{m+i}+\dfrac{1}{2}a_{i}A_{0}+\dfrac{1}{2}\sum_{m=1}^{i-1}a_{m}A_{i-m}+\dfrac{1}{2}\sum_{m=i+1}^{h}a_{m}A_{m-i}+$
$\displaystyle+\dfrac{1}{2}\sum_{m=1}^{h}b_{m}B_{m+i}+\dfrac{1}{2}b_{i}B_{0}-\dfrac{1}{2}\sum_{m=1}^{i-1}b_{m}B_{i-m}+\dfrac{1}{2}\sum_{m=i+1}^{h}b_{m}B_{m-i}=$
$\displaystyle=a_{0}A_{i}+a_{i}A_{0}+\dfrac{1}{2}\sum_{m=1}^{h-i}\left(a_{m}A_{m+i}+b_{m}B_{m+i}\right)+\dfrac{1}{2}\sum_{m=1}^{i-1}\left(a_{m}A_{i-m}-b_{m}B_{i-m}\right)+$
$\displaystyle+\dfrac{1}{2}\sum_{m=i+1}^{h}\left(a_{m}A_{m-i}+b_{m}B_{m-i}\right),$
$\displaystyle\beta_{i}$
$\displaystyle=a_{0}B_{i}+\dfrac{1}{2}\sum_{m=1}^{\infty}a_{m}B_{m+i}-\dfrac{1}{2}\sum_{m=1}^{\infty}a_{m}B_{m-i}-\dfrac{1}{2}\sum_{m=1}^{\infty}b_{m}A_{m+i}+\dfrac{1}{2}\sum_{m=1}^{\infty}b_{m}A_{m-i}=$
$\displaystyle=a_{0}B_{i}+\dfrac{1}{2}\sum_{m=1}^{h-i}a_{m}B_{m+i}+\dfrac{1}{2}\sum_{m=1}^{i-1}a_{m}B_{i-m}-\dfrac{1}{2}\sum_{m=i+1}^{h}a_{m}B_{m-i}-$
$\displaystyle-\dfrac{1}{2}\sum_{m=1}^{h-i}b_{m}A_{m+i}+b_{i}A_{0}+\dfrac{1}{2}\sum_{m=1}^{i-1}b_{m}A_{i-m}+\dfrac{1}{2}\sum_{m=i+1}^{h}b_{m}A_{m-i}=$
$\displaystyle=a_{0}B_{i}+b_{i}A_{0}+\dfrac{1}{2}\sum_{m=1}^{h-i}\left(a_{m}B_{m+i}-b_{m}A_{m+i}\right)+\dfrac{1}{2}\sum_{m=1}^{i-1}\left(a_{m}B_{i-m}+b_{m}A_{i-m}\right)+$
$\displaystyle+\dfrac{1}{2}\sum_{m=i+1}^{h}\left(-a_{m}B_{m-i}+b_{m}A_{m-i}\right).$

Applying the obtained formulas to calculate the products of trigonometric
polynomials to the residuals, we can write the equations for the $i$-th
harmonics ($i=\overline{1,h}$ is the number of harmonics, $k=\overline{1,3}$
is residual number):

$k=1$:

$\begin{array}[]{r}i\omega s_{1,i}-10c_{2,i}+10c_{1,i}=0,\vskip 6.0pt plus
2.0pt minus 2.0pt\\\ -i\omega c_{1,i}-10s_{2,i}+10s_{1,i}=0,\end{array}$

the equation corresponding to the constant term for the first residual is

$x_{1,0}-x_{2,0}=0,$

$k=2$:

$\displaystyle i\omega
s_{2,i}-28c_{1,i}+c_{2,i}+x_{1,0}c_{3,i}+c_{1,i}x_{3,0}$
$\displaystyle+\dfrac{1}{2}\sum_{m=1}^{h-i}\left(c_{1,m}c_{3,m+i}+s_{1,m}s_{3,m+i}\right)+$
$\displaystyle+\dfrac{1}{2}\sum_{m=1}^{i-1}\left(c_{1,m}c_{3,i-m}-s_{1,m}s_{3,i-m}\right)+$
$\displaystyle+\dfrac{1}{2}\sum_{m=i+1}^{h}\left(c_{1,m}c_{3,m-i}+s_{1,m}s_{3,m-i}\right)=0,$
$\displaystyle-i\omega
c_{2,i}-28s_{1,i}+s_{2,i}+x_{1,0}s_{3,i}+s_{1,i}x_{3,0}$
$\displaystyle+\dfrac{1}{2}\sum_{m=1}^{h-i}\left(c_{1,m}s_{3,m+i}-s_{1,m}c_{3,m+i}\right)+$
$\displaystyle+\dfrac{1}{2}\sum_{m=1}^{i-1}\left(c_{1,m}s_{3,i-m}+s_{1,m}c_{3,i-m}\right)+$
$\displaystyle+\dfrac{1}{2}\sum_{m=i+1}^{h}\left(-c_{1,m}s_{3,m-i}+s_{1,m}c_{3,m-i}\right)=0,$

the equation corresponding to the constant term for the second residual is

$-28x_{1,0}+x_{2,0}+x_{1,0}x_{3,0}+\dfrac{1}{2}\sum_{m=1}^{h}\left(c_{1,m}c_{3,m}+s_{1,m}s_{3,m}\right)=0,$

$k=3$:

$\displaystyle i\omega s_{3,i}-x_{1,0}c_{2,i}-c_{1,i}x_{2,0}$
$\displaystyle-\dfrac{1}{2}\sum_{m=1}^{h-i}\left(c_{1,m}c_{2,m+i}+s_{1,m}s_{2,m+i}\right)-$
$\displaystyle-\dfrac{1}{2}\sum_{m=1}^{i-1}\left(c_{1,m}c_{2,i-m}-s_{1,m}s_{2,i-m}\right)-$
$\displaystyle-\dfrac{1}{2}\sum_{m=i+1}^{h}\left(c_{1,m}c_{2,m-i}+s_{1,m}s_{2,m-i}\right)+\dfrac{8}{3}c_{3,i}=0,$
$\displaystyle-i\omega c_{3,i}-x_{1,0}s_{2,i}-s_{1,i}x_{2,0}$
$\displaystyle-\dfrac{1}{2}\sum_{m=1}^{h-i}\left(c_{1,m}s_{2,m+i}-s_{1,m}c_{2,m+i}\right)-$
$\displaystyle-\dfrac{1}{2}\sum_{m=1}^{i-1}\left(c_{1,m}s_{2,i-m}+s_{1,m}c_{2,i-m}\right)-$
$\displaystyle-\dfrac{1}{2}\sum_{m=i+1}^{h}\left(-c_{1,m}s_{2,m-i}+s_{1,m}c_{2,m-i}\right)+\dfrac{8}{3}s_{3,i}=0,$

the equation corresponding to the constant term for the third residual is

$-x_{1,0}x_{2,0}-\dfrac{1}{2}\sum_{m=1}^{h}\left(c_{1,m}c_{2,m}+s_{1,m}s_{2,m}\right)+\dfrac{8}{3}x_{3,0}=0,$

the additional system equation is

$x_{3,0}+\sum_{i=1}^{h}c_{3,i}-27=0.$

##  5 The Results of the Computational Experiment

As a result of numerous computational experiments, the initial approximation
was chosen for the cyclic frequency, constant terms, and amplitudes at
$h=h_{1}=5$:

$\begin{array}[]{c}\omega=4,\>\>x_{1,0}=x_{2,0}=x_{3,0}=0,\>\>c_{1,i}=-1,\>i=\overline{1,5},\\\
s_{1,j}=0,\>j=1,3,4,5,\>\>s_{1,2}=1.\end{array}$

This result is remarkable in that the Newton method converges to a solution
different from the equilibrium positions. Therefore, to improve the accuracy
of the approximate periodic solution, we consider a system of algebraic
equations for the value of $h$ equal to some $h_{2}>h_{1}$.

The obtained numerical solution of the system for $h=h_{1}$ is taken as the
initial approximation for amplitudes with indices $i\leq h_{1}$ for a system
with $h=h_{2}$, and the values of the initial approximation for amplitudes
with indices $i>h_{1}$ are assumed to be zero.

Table 1: The amplitudes of harmonics for $\tilde{x}_{1}(t)$, $x_{1,0}=0$. $i$ | $c_{1,i}$ | $s_{1,i}$
---|---|---
1 | $-5.780478259196228$ | $8.56017654325353$
2 | 0 | 0
3 | $3.160762628380509$ | $2.239212141102876$
4 | 0 | 0
5 | $0.6958870387616096$ | $-0.7979388979225431$
6 | 0 | 0
7 | $-0.1891992374027477$ | $-0.1864921358925765$
8 | 0 | 0
9 | $-0.04770429623010056$ | $0.04554044367245914$
10 | 0 | 0
11 | $0.01112322884679491$ | $0.01209138588669679$
12 | 0 | 0
13 | $0.003061207095371694$ | $-0.002735092350544739$
14 | 0 | 0
15 | $-6.744578887916229\cdot 10^{-4}$ | $-7.748319471034087\cdot 10^{-4}$
16 | 0 | 0
17 | $-1.960718247379475\cdot 10^{-4}$ | $1.665584161919807\cdot 10^{-4}$
18 | 0 | 0
19 | $4.116738805347028\cdot 10^{-5}$ | $4.960493476144467\cdot 10^{-5}$
20 | 0 | 0
21 | $1.254757391175977\cdot 10^{-5}$ | $-1.018054283421179\cdot 10^{-5}$
22 | 0 | 0
23 | $-2.518375902000733\cdot 10^{-6}$ | $-3.173486439630506\cdot 10^{-6}$
24 | 0 | 0
25 | $-8.025338211960923\cdot 10^{-7}$ | $6.230623750431923\cdot 10^{-7}$
26 | 0 | 0
27 | $1.541534734542893\cdot 10^{-7}$ | $2.0292802821633\cdot 10^{-7}$
28 | 0 | 0
29 | $5.130649139299358\cdot 10^{-8}$ | $-3.813725452268523\cdot 10^{-8}$
30 | 0 | 0
31 | $-9.43393531993558\cdot 10^{-9}$ | $-1.297038481588497\cdot 10^{-8}$
32 | 0 | 0
33 | $-3.278552746800046\cdot 10^{-9}$ | $2.333260259021725\cdot 10^{-9}$
34 | 0 | 0
35 | $5.76957885768651\cdot 10^{-10}$ | $8.28626640138045\cdot 10^{-10}$
Table 2: The amplitudes of harmonics for $\tilde{x}_{2}(t)$, $x_{2,0}=0$. $i$ | $c_{2,i}$ | $s_{2,i}$
---|---|---
1 | $-2.32972926505593$ | $10.89038310357172$
2 | 0 | 0
3 | $5.86875317198698$ | $-1.5832552129833$
4 | 0 | 0
5 | $-0.9124249133801483$ | $-2.200556873678218$
6 | 0 | 0
7 | $-0.7154457265566421$ | $0.3473932955614448$
8 | 0 | 0
9 | $0.1175186702136983$ | $0.2186139734768588$
10 | 0 | 0
11 | $0.06473984670858603$ | $-0.03723215039412078$
12 | 0 | 0
13 | $-0.01127208646321726$ | $-0.01877739524860192$
14 | 0 | 0
15 | $-0.005359671824365359$ | $0.003303445299126894$
16 | 0 | 0
17 | $9.453499475830811\cdot 10^{-4}$ | $0.001510235036151227$
18 | 0 | 0
19 | $4.211022386354685\cdot 10^{-4}$ | $-2.657049331814368\cdot 10^{-4}$
20 | 0 | 0
21 | $-7.363528144366622\cdot 10^{-5}$ | $-1.164013765469982\cdot 10^{-4}$
22 | 0 | 0
23 | $-3.19419300699788\cdot 10^{-5}$ | $2.017609175377016\cdot 10^{-5}$
24 | 0 | 0
25 | $5.47663534401654\cdot 10^{-6}$ | $8.710929378319451\cdot 10^{-6}$
26 | 0 | 0
27 | $2.362852034076972\cdot 10^{-6}$ | $-1.474901091428546\cdot 10^{-6}$
28 | 0 | 0
29 | $-3.94532524722541\cdot 10^{-7}$ | $-6.379296603810031\cdot 10^{-7}$
30 | 0 | 0
31 | $-1.715198229248314\cdot 10^{-7}$ | $1.049218598356554\cdot 10^{-7}$
32 | 0 | 0
33 | $2.776045093375681\cdot 10^{-8}$ | $4.59473450493284\cdot 10^{-8}$
34 | 0 | 0
35 | $1.22681173575872\cdot 10^{-8}$ | $-7.31171826830086\cdot 10^{-9}$
Table 3: The amplitudes of harmonics for $\tilde{x}_{3}(t)$, $x_{3,0}=23.04210397942006$. $i$ | $c_{3,i}$ | $s_{3,i}$
---|---|---
1 | 0 | 0
2 | $7.568410271550653$ | $-9.50386584559212$
3 | 0 | 0
4 | $-3.555327211552558$ | $-1.844710563805469$
5 | 0 | 0
6 | $-0.4741220131932616$ | $1.279043179069961$
7 | 0 | 0
8 | $0.4227292179138024$ | $0.1274574086305204$
9 | 0 | 0
10 | $0.03498415351761577$ | $-0.1315337800809524$
11 | 0 | 0
12 | $-0.03934013541135439$ | $-0.009645786231708874$
13 | 0 | 0
14 | $-0.002660052258813564$ | $0.01145537653603837$
15 | 0 | 0
16 | $0.003271688724557337$ | $7.33752523103949\cdot 10^{-4}$
17 | 0 | 0
18 | $2.024982256871223\cdot 10^{-4}$ | $-9.206266886554897\cdot 10^{-4}$
19 | 0 | 0
20 | $-2.560063570343799\cdot 10^{-4}$ | $-5.58964460662525\cdot 10^{-5}$
21 | 0 | 0
22 | $-1.542436654918173\cdot 10^{-5}$ | $7.050327849098175\cdot 10^{-5}$
23 | 0 | 0
24 | $1.926014222030195\cdot 10^{-5}$ | $4.25261452471065\cdot 10^{-6}$
25 | 0 | 0
26 | $1.170939944189529\cdot 10^{-6}$ | $-5.225643926851625\cdot 10^{-6}$
27 | 0 | 0
28 | $-1.409525591131397\cdot 10^{-6}$ | $-3.21879984959824\cdot 10^{-7}$
29 | 0 | 0
30 | $-8.83134288999026\cdot 10^{-8}$ | $3.782652721710986\cdot 10^{-7}$
31 | 0 | 0
32 | $1.010610960272394\cdot 10^{-7}$ | $2.418021923473667\cdot 10^{-8}$
33 | 0 | 0
34 | $6.606163280924149\cdot 10^{-9}$ | $-2.689431432873997\cdot 10^{-8}$
35 | 0 | 0
Figure 1: The cycle obtained by described method.

Tables 1–3 show the result of solving the system for $h=35$; the accuracy of
the Newton method is $10^{-8}$. The period value is obtained equal to
$T=1.558652210$, the initial condition for the obtained approximate periodic
solution is

$\begin{array}[]{c}\tilde{x}_{1}(0)=-2.147367631,\>\>\tilde{x}_{2}(0)=2.078048211,\>\>\tilde{x}_{3}(0)=27.\end{array}$
(6)

The initial values (6) were checked on the period in a computer program that
implements the numerical integration of the system (1) by the modified power
series method [8] with an accuracy of estimating the common term of the series
$10^{-25}$, 100 bits for mantissa real number and machine epsilon
$1.57772\cdot 10^{-30}$. With such parameters of the method, the approximate
values of the phase coordinates obtained by numerical integration were also
verified by the same numerical method, but in reverse time. The values in the
reverse time coincide with (6) up to the 9th character inclusive after the
point. The resulting values of $x_{1}(T)$, $x_{2}(T)$ and $x_{3}(T)$ coincide
with (6) up to the 8th character inclusive.

The cycle corresponding to (6) is shown in Fig. 1. Note that the cycle found
coincides with the first Viswanath cycle in [6], all signs after the point for
$T$ also coincide with the data from [6].

##  6 Acknowledgements

The reported study was funded by RFBR according to the research project
20-01-00347.

## References

  * [1] Lorenz, E. N. Deterministic Nonperiodic Flow, Journal of the Atmospheric Sciences, vol. 20, no. 2 (1963), pp. 130-141.
  * [2] Tucker, W. A Rigorous ODE Solver and Smale’s 14th Problem, Foundations of Computational Mathematics, vol. 2, no. 1 (2002), pp. 53-117.
  * [3] Rabinovich, M. I. Stochastic Self-Oscillations and Turbulence, Soviet Physics Uspekhi, vol. 21, no. 5 (1978), pp. 443-469.
  * [4] Galias, Z., Tucker, W. Validated Study of the Existence of Short Cycles for Chaotic Systems Using Symbolic Dynamics and Interval Tools, International Journal of Bifurcation and Chaos, vol. 21, no. 2 (2011), pp. 551-563.
  * [5] Lozi, R. Can We Trust in Numerical Computations of Chaotic Solutions of Dynamical Systems?, Topology and Dynamics of Chaos. In Celebration of Robert Gilmore’s 70th Birthday. - World Scientific Series in Nonlinear Science Series A, vol. 84 (2013), pp. 63-98.
  * [6] Viswanath, D. The Fractal Property of the Lorenz Attractor, Physica D: Nonlinear Phenomena, vol. 190, no. 1-2 (2004), pp. 115-128.
  * [7] Viswanath, D. The Lindstedt-Poincare Technique as an Algorithm for Computing Periodic Orbits, SIAM Review, vol. 43, no. 3 (2001), pp. 478-495.
  * [8] Pchelintsev, A. N. Numerical and Physical Modeling of the Dynamics of the Lorenz System, Numerical Analysis and Applications, vol. 7, no. 2 (2014), pp. 159-167.
  * [9] Neymeyr, K., Seelig, F. Determination of Unstable Limit Cycles in Chaotic Systems by Method of Unrestricted Harmonic Balance, Zeitschrift für Naturforschung A, vol. 46, no. 6 (1991), pp. 499-502.
  * [10] Luo, A. C. J., Huang, J. Approximate Solutions of Periodic Motions in Nonlinear Systems via a Generalized Harmonic Balance, Journal of Vibration and Control, vol. 18, no. 11 (2011), pp. 1661-1674.
  * [11] Luo, A. C. J. Toward Analytical Chaos in Nonlinear Systems, John Wiley & Sons, Chichester, ISBN: 978-1-118-65861-1, 2014, 258 pp.
  * [12] Luo, A. C. J., Guo, S. Analytical Solutions of Period-1 to Period-2 Motions in a Periodically Diffused Brusselator, Journal of Computational and Nonlinear Dynamics, vol. 13, no. 9, 090912 (2018), 8 pp.
  * [13] Pchelintsev, A. N. The Programs for Finding of Periodic Solutions in the Lorenz Attractor, GitHub, https://github.com/alpchelintsev/periodic_sols
  * [14] Tolstov, G. P. Fourier Series, Dover Publications, New York (1962), 336 pp.