Principal Component Analysis in the Nonlinear Dynamics of Beams : Purification of the Signal from Noise Induced by the Nonlinearity of Beam Vibrations

1Department of Applied Mathematics and Systems Analysis, Saratov State Technical University, 77 Politeknicheskaya Str., Saratov 41054, Russia 2Cybernetics Institute, National Research Tomsk Polytechnic University, 30 Lenin Avenue, Tomsk 634050, Russia 3Department of Automation, Biomechanics and Mechatronics, Lodz University of Technology, 1/15 Stefanowski St., 90-924 Lodz, Poland 4Department of Mathematics and Modelling, Saratov State Technical University, 77 Politeknicheskaya Str., Saratov 41054, Russia


Introduction
Numerous methods can be employed to clean a signal from noise.For instance, wavelet transform [1], Fourier transform [2], Hilbert-Huang transform [3], or empirical mode decomposition [4] is often applied for this purpose.In all of the methods, the goal is to decompose the signal to a variety of fundamental constituents and enable analysing each of them separately.As a result, various features of the signal can be detected and the signal can be properly analysed.In this paper, the attention is paid to the method called the principal component analysis (PCA), which is widely used to amplify variation in the spectral data and reveal strong patterns in a dataset.
PCA [5] is an approach usually utilised to simplify investigation and visualisation of data.In particular, it is commonly used in pattern recognition and compression of various types of data, especially images, having high correlation between their components.For instance, it is implemented in frontal face databases, where the goal is to process images in such a way that noise is removed from a neighbourhood of a block of pixels.In order to accomplish it, it is necessary to present this block as a set of points in a multidimensional space.Then, the PCA can be applied, as a result of which only the first conversion components, assumed to contain the most useful information, are left.The remaining components are said to comprise unwanted noise.Thus, if an inverse transformation is applied to the first conversion components, a denoised image is obtained as output.One of advantages of this method is that no specified accuracy is required.
To assess the number of principal components in the appropriate proportion of variance, an objective approach can be always applied.However, when the signal and noise are not separated, no predetermined accuracy is of particular 2 Advances in Mathematical Physics 2h z a q(x, t) = q 0 sin( p t) x meaning.Thus, if one considers projections of the components on different planes, the components of the signal are tightly packed, keeping the amplitude relatively large.On the contrary, noise components associated with a relatively small amplitude occupy larger space (are scattered more strongly).Thus, in the above-described case, the principal component analysis plays a role of a filter.Namely, the signal is mainly contained in the projection of the first principal components, whereas other components exhibit much higher noise content.
In this paper, the principal component analysis is employed to analyse the signal coming from vibrations of beams.For this purpose, it is necessary to supplement this method with a mathematical model of a vibrating beam.As it is well known, there are several theories associated with the beam model.Particularly, Bernoulli-Euler [6] and Timoshenko model [7] are widely employed in many mechanical problems.In addition, a generalisation of the Timoshenko model, invented in 1964 by two Ukrainian scientists, Sheremetev and Pelekh [8], is often considered.It should be mentioned that this model was rediscovered 27 years later by Levinson [9] and Reddy [10]; thus in English literature, it is typically called the Levinson-Reddy model.
In this paper, a comparative analysis of the amplitudefrequency characteristics for all the above-mentioned, that is, Euler-Bernoulli, Timoshenko, and Sheremetev-Pelekh-Levinson-Reddy, beam models is conducted.Comparison of the amplitude-frequency characteristics is based on the signal analysis of solutions of systems of linear and nonlinear differential equations.

Problem Formulation
The study is conducted on a single-layer beam occupying a two-dimensional region of  2 space with the Cartesian coordinate system  introduced in the following way: the  axis is directed from the left to the right along the beam midline, and the  axis is directed downwards, perpendicularly to the  axis (see Figure 1).In the introduced coordinate system, the two-dimensional domain Ω of the beam is defined as follows: Ω = { ∈ [0, ]; −ℎ ≤  ≤ ℎ}, 0 ≤  ≤ ∞.Here and below we use the following notation: ℎ denotes the beam height and  stands for the beam length.The beam is subjected to the load (, ) acting in the direction normal to the midline of the beam.In addition, the beam lies on an elastic Winkler foundation.
The mathematical model of the Sheremetev-Pelekh-Levinson-Reddy beam is a third-order approach based on the following hypotheses: (i) rotational inertia effects of beam elements are ignored, but inertial forces responsible for the displacement of the beam along the line normal to the midline are taken into consideration; (ii) external forces do not change their direction during deformations of the beam; (iii) the beam longitudinal dimension is significantly larger than its transverse dimensions; (iv) geometric nonlinearity is taken in the form proposed by von Kármán; (v) a beam cross-section does not remain flat and perpendicular to the deformed beam axis; (vi) the inertial components associated with rotation are included.
The displacements of an arbitrary beam point are as follows: where   is the angle of rotation of the normal to the beam midline,   and   are unknown functions, and  stands for the beam normal deflection.The equation of the beam motion and boundary conditions for this equation can be obtained using a standard approach, that is, by computing kinetic and potential energies as well as works done by external forces [11].
The true trajectory differs from other possible trajectories when the following condition is satisfied: where  and ^are variations of the kinetic and potential energies of the system, respectively, and    is the sum of elementary works of external forces.Using the variational principles, we obtain the equations of motion taking into account the energy dissipation.
A mathematical model of the beam based on the Sheremetev-Pelekh-Levinson-Reddy hypothesis and including the geometric nonlinearity is governed by a system of the following nonlinear PDEs: where , and  4 (, ) = (/)( 2 / 2 ) are nonlinear operators, (, ) is the beam deflection in the normal direction, (, ) is the beam displacement in the beam longitudinal direction,  is damping coefficient, E is Young's modulus, ℎ is beam height,  is specific material gravity,  is the gravity of Earth,  is time, and  1 is a proportionality constant representing contact pressure per unit beam deformation, commonly referred to as the modulus of subgrade reaction [12,13].Thus,  is the only quantity characterising the functionally graded material.
The below dimensionless parameters are introduced: Also, the following boundary condition (fixed support) is added to (3): and the following initial conditions are taken: The above equations contain a fourth-order derivative with respect to a coordinate variable, which is extremely important in proving the existence of solutions to the mentioned PDEs as well as the convergence of different methods employed to find them.
The Sheremetev-Pelekh-Levinson-Reddy beam is governed by a system of linear PDEs: Boundary and initial conditions ( 5) and ( 6) are also employed, respectively.
The mathematical model of the Timoshenko beam assumes that the normal does not remain perpendicular to the beam midline and it rotates by an angle   , remaining rectilinear.The formulas for the displacements of an arbitrary beam point have the form In this case, the following system of governing PDEs is derived: where the dimensionless parameters introduced in formulas (4) have been already employed.Equations ( 9) are supplemented with the following boundary condition (fixed support): and initial conditions (6).
The linear PDEs governing the Timoshenko beam follow ODEs (11) are supplemented with the boundary ( 9) and initial (6) conditions.Note that linear Timoshenko PDEs are of the second order with respect to , which sometimes makes proving the convergence of the employed numerical methods difficult.
On the other hand, the Euler-Bernoulli hypothesis assumes that the line perpendicular to the beam midline before the deformation process remains perpendicular also during deformation but rotates by the angle   .Using formula (8) and replacing the shift function   with the angle of rotation − 2 / 2 , the following Euler-Bernoulli equation is obtained: Equation ( 11) is already in the dimensionless form, in which a part of relations ( 4) is substituted with For simplicity, bars over dimensionless parameters have been already omitted in (11).
In the case of the fixed support, we supplement the equations with the following: (i) boundary conditions: (ii) initial conditions: The counterpart linear PDE takes the following form: and the boundary conditions ( 14) and the initial conditions ( 15) are employed here.

Solution Methods
The same solving algorithm can be used for all the sofar defined problems.It is clear that (3)-( 16) cannot be solved analytically, and hence they are solved numerically by means of applying the finite-difference method (FDM) of the second-order accuracy with respect to the spatial coordinate  [14].In order to reduce PDEs to ODEs with respect to time, we use the finite-difference approximations applying Taylor expansion in the neighbourhood of a point   .Consider a mesh The following difference operators including the approximation of ( 2 ), where  is the step of spatial coordinate partition, are introduced: Hence, PDEs are reduced to the second-order ODEs with respect to time with corresponding boundary and initial conditions, based on the application of the finitedifference approximation of the second-order accuracy.Then, the second-order ODEs are reduced to a first-order system, which is solved by the fourth-order Runge-Kutta method.
For example, a system of nonlinear PDEs (3), ( 5), and ( 6) is reduced to the following ODEs: with the associated boundary conditions Validity of the results for all considered problems is illustrated and discussed in [7].

The Method of Principal Components.
The principal component analysis is a procedure of variable reduction.It is particularly useful when one analyses data of a (possibly large) number of variables and believes there is some redundancy in those variables; that is, some of the variables are correlated with one another, possibly as a result of measuring the same feature.Since a part of the data is redundant, it is believed that it should be possible to reduce the observed variables into a smaller number of principal components (artificial variables) that will account for the most of the variance in the observed variables.In other words, the physical meaning of this method is that the multidimensional data matrix can be reduced so that a number of possibly correlated variables are transformed into a smaller (or the same) number of linearly uncorrelated, orthogonal variables.After applying the PCA, the variance of the transformed data is maximised.
The method of principal component analysis cleanses the signal from noise and localizes the main frequency signal.The input data matrix is organised in such a way the variables are placed in columns while samples are in rows.To separate the given signal from noise, the linear matrix decomposition   =      +   is employed, where the initial   = (  ) , =1,=1 data matrix is purified from the noise. denotes the number of principal components,   = (  ) ,= =1,=1 is called the score matrix,   = (  ) ,= =1,=1 is called the loading matrix, and   = (  ) , =1,=1 stands for the residual matrix [15].To conduct the decomposition, first of all, the matrix is centred, which means that mean values of each column of the matrix are subtracted from these columns, respectively.
Then, a covariance matrix     T is usually calculated, and the eigenvalues and eigenvectors are found.After arranging the eigenvalues in the high-to-low manner, the eigenvalues should be projected onto a new space.The eigenvector of the highest eigenvalue is said to be the first principal component, the one of the second highest eigenvalue, as the second principal component, and so forth.
To find the eigenvalues and eigenvectors, also the singular-value decomposition (SVD) is often employed.Its detailed description can be found in [16,17], where also the detailed descriptions of the mathematical and physical meaning of spurious eigenvalues and fictitious frequencies in a nontrivial boundary solution are studied based on a SVD updating technique.
In SVD, a real matrix can be decomposed to the product of the so-called left unitary matrix , the diagonal matrix  with positive elements called the singular values of   , and the transpose of the right unitary matrix   .The SVD results in the eigenvalues and eigenvectors needed for PCA (the squared eigenvalues are equal to those found by covariance matrix).
Also, a so-called truncated SVD method can be distinguished, in which only a chosen number (corresponding to the largest singular values of ) of vectors of  and   are studied.Hence, the number of computational operations can be reduced.

Selection of the Number of Principal Components.
In what follows, a selection of the number of the principal components is considered.First, the resulting matrix is centred.From each column is subtracted, where  stands for a number of rows.Then, the matrix is standardised, that is, each column   is divided by its standard deviation: The combination of centring and normalisation by columns is called autoscaling, which yields Note that average values of many variables differ from zero.In addition, also standard deviations are substantially different.However, after the autoscaling procedure, the average of all variables becomes equal to zero and the deviation becomes equal to one.Let us represent the matrix W = (w  ) , =1,=1 employing a linear decomposition W =   + .To calculate the score matrix  and loadings matrix , the singular-value decomposition (SVD) can be used.The SVD of the matrix W obtained after autoscaling yields W =   , where  =  ( = (  ) , =1,=1 ),  =   .Note that matrices  = (  ) , =1,=1 and  = (V  ) , =1,=1 are orthogonal and that the main diagonal of the matrix  contains the eigenvalues, whereas its other elements are equal to zero.
The value of the principal component  is chosen so that the eigenvalues of the matrix  are greater than one or that the dependences of the total residual variance (TRV) and the explained residual variance (ERV) on the number of principal components , that is, TRV() and ERV(), dramatically change their behaviour.Calculation of TRV() and ERV() is described in the below sections.Once the number of principal components is chosen, then the residual matrix   = (  ) , =1,=1 is yielded by the formula:   = W −      .The matrix has a variable dimension   = (  ) , =1,=1 and   = (  ) , =1,=1 .The important property of the principal component analysis (PCA) is the orthogonality (independence) of principal components.
Therefore, the score matrix   is not reconstructed by increasing the number of components, but another column corresponding to the new direction is just added.The same holds with the matrix of loadings   .
To find TRV and ERV, it is necessary to compute the residual matrix:   = (  ) , =1,=1 for each number of principal components  as well as to compute the required indicators.
The following computational steps are required: (1) Find ]  = ∑ = =1    , which determines the square of the distance between the original vector   and its projection on the PC space.The smaller the value is, the better the vector approaches the th sample.
(2) Calculate (for all samples) the average square of the distance as

Numerical Results
Let us consider free vibrations of a fixed beam on the elastic Winkler foundations with the subgrade modulus  1 = 1.The beam is subjected to uniformly distributed transverse load  =  0 sin(  ) with the amplitude  0 = 100 and the frequency   = 5.The dissipation coefficient  = 0.The ratio of the beam length  to its thickness 2ℎ is  = /2ℎ = 30.Table 1 shows graphs of power spectra constructed with the help of Fast Fourier Transform.
4.1.The Influence of Geometric Nonlinearity on the Amplitude-Frequency Characteristics of the Signal.Let us investigate the effect of geometric nonlinearity on the solution of the studied PDEs.A signal of a small amplitude within the framework of linear theory for thin beams ( ≤ 0.25(2ℎ)) and the influence of geometric nonlinearity on the frequency spectrum of the signal at small deflection of the beam have been studied.
By solving linear PDEs (11), (10), and (6) taking into account Timoshenko hypothesis, it has been detected that quasi-periodic vibrations of the beam occur at four frequencies, that is, the excitation frequency   = 5 and linearly dependent frequencies:  1 = (1/5)  ;  2 = 2 1 ;  3 =   −  1 .By taking into account the geometric nonlinearity of the beam midline deformation (8), (9), and (5), one can follow the qualitative changes of the beam vibrations: In what follows, the amplitude-frequency characteristics of the signal obtained by solving a system of linear PDEs are considered (Table 1, the first column).The results show that the beam vibrates at four frequencies: the frequency of the excitation   = 5 and linearly dependent frequencies  1 = (1/5)  ;  2 = 2 1 ;  3 =   −  1 .In linear systems, solutions to linear PDEs based on the Timoshenko and Sheremetev-Pelekh-Levinson-Reddy hypotheses (7), (5), and (6) give the same frequency characteristics for all the studied signals.
Parameters (features) TRV and ERV indicate what percentage of noise remains after projection of the signal on the multidimensional space PC1 − PCA.In other words, these characteristics indicate what number of principal components is necessary to cleanse the signal from noise.
) Compute the total residual variance TRV = V 0 / and the explained residual variance ERV = 1 − ( ⋅ V 0 )/ ∑  w 2  .In order to separate the signal from noise, the linear matrix decomposition   =      +   is used, where the matrix   = (  ) , =1,=1 is separated from noise and where  stands for the number of principal components.

Table 1 :
Power spectra of signals obtained for linear and nonlinear cases, constructed for different beam models.