1. Introduction
Acoustic signal processing is a distinctive application in numerous sectors in life, like industry and communications. In electric motor diagnosis, acoustic noise investigation is a robust alternative to complement the information given by other methods [1,2,3,4]. A common issue in these machines is the occasional occurrence of rotor damage, such as broken bars. These faults may indeed lead to disastrous disappointments and constrained motor blackouts that can infer critical misfortunes for the businesses using them. In particular, it is especially important in the case of high-voltage motors that are used in exceptionally large machines. Figure 1 depicts a picture of the rotor of a motor with broken bars.
The acoustic noise analysis permits rotor damage to be detected and adjacent broken rotor bars to be classified, as proven in previous works [5]. With regards to the detection of non-adjacent bar breakages, it has not been reliably solved by other techniques like vibration or current analysis techniques that have even provided false-negative diagnostics when detecting this fault [6].
Noise-based diagnosis algorithms based on fast Fourier transform (FFT) and other spectral methods that rely on subspace vectors were considered in [5,6,7,8,9,10,11,12]. However, the FFT has a drawback: it is not immune to noise since the spectrum of a noisy signal also includes the noise spectrum. To overcome this constraint, a few researchers have utilized high-resolution spectral algorithms, with the impediment of not knowing a priori the number of subspaces designated to the noise [7,8].
Recently, a filter based on the windowed Fourier transforms (WFT) [13] has also been considered, a strategy that lies in a choice of amplitudes: multiexpanded in the spectrum [14], or an algorithm based on calculations on spectral vector subspaces such as MUSIC and ESPRIT [15]. Furthermore, wavelet-based strategies [11,16,17,18,19,20,21,22] or empirical mode decomposition (EMD) [23] have also been considered.
The application of higher-order cumulants or bispectrum, or its one-dimensional component, has been employed for fault diagnosis in [24,25,26,27]. However, the use of all the data within the bispectrum lies in two-dimensional Fourier transforms that may increase computational cost.
Other methods different from the examination of sound signals allow the number of broken bars to be distinguished and the deficiencies in a motor to be analyzed, such as artificial intelligence algorithms [12,28,29,30], thermal images processing [6], and chaos theory methods [31].
Nevertheless, none of them have demonstrated being substantially sufficient for the case of non-adjacent broken bar diagnosis. Subsequently, we propose a method based on calculations that imply second-order statistics, convolution, and spectral subtractions for the failure detection of bar breakages, using the noise of an induction motor. These calculations are used not only as it were to identify adjacent bar breakages, but also to examine the plausibility of identifying non-adjacent broken bars and hence improve the outcomes of classical methods.
As commented, one of the cornerstones of the condition monitoring area is the detection of non-adjacent broken bars. Diagnosis tools like the conventional motor current signature analysis (MCSA) may lead to false-positive detection of non-adjacent breakages [32]. Recent efforts have been addressed to detect rotor bar failures, regardless of if the bars were consecutive or not [33,34]. Here, we propose a new acoustic approach that enables us to detect failures. The method also works for the detection of non-adjacent broken bars.
In Section 2, we present the foundations of the analyzed fault and of the tools considered in our work. In Section 3, we show that the convolution of a signal motor with its autocovariance permits diminishing the noise. We present the proposed design and the results in Section 4, and the conclusions in Section 5.
2. Materials and Methods
2.1. Faults Analyzed: Broken Rotor Bars (BRB)
It is already known that a broken bar in a motor leads to a mechanical distortion in the air gap magnetic field. This distortion yields certain harmonics in the stator phase current. The foremost one is the lower sideband harmonic (also known as left sideband harmonic), with frequency given by , where stands for the line supply frequency, for the slip defined in (1), for the synchronous speed of the machine (in r.p.m., equivalent to , with numbers of poles), and for the motor speed.
(1)
We know that the lower sideband harmonic is linked to a torque oscillation, and it yields another harmonic that appears in the stator current spectrum: the upper sideband harmonic (or right sideband harmonic), given by
Moreover, the speed oscillation caused by the presence of broken rotor bars generates a frequency modulation in the rotational frequency, yielding two sidebands in the vibration spectrum and, hence, in the noise spectrum, with frequencies given by (2) [23].
(2)
2.2. The Noise-Reduction Algorithm
Concerning the processing of signals for noise reduction, different works have been carried out. These works can be classified according to various criteria. However, according to the number of signals at the input of the noise-reduction system, they can be classified into methods that use two or more input signals, and methods that use only one input signal [22]. A classification proposal can then be provided; it tries to group different techniques according to their characteristics and requirements:
Classical FIR (Finite Impulse Response) filtering methods [35,36,37].
Adaptive noise reduction methods [38,39,40].
Artificial intelligence methods. Neural Networks [41,42].
Wavelet-based methods [43,44,45].
Statistical signal processing methods [46,47].
From previous works, classical filtering techniques limit their use to reduce noise only in the frequency band to which the filter is limited. On the other hand, although adaptive filtering algorithms have the ability to adapt, most of them need a reference sample for obtaining tangible results. Similarly, artificial intelligence and wavelet-based methods need a reference to adjust the desired output and a comparison threshold, which depends on the noise power, respectively. The proposed noise reduction algorithm is intended to simplify the signal and keep only the information relative to the motor behavior that leads to identifying the broken bars. It distinguishes and separates the spectral components of the acoustics that do not inform about them. These components are considered as random interferences from the environment. We consider them part of a Gaussian noise. We split the propose method into the following steps:
We set a filter from the convolution of the signal with its autocovariance. More details about these operations will be shown below and in the recent survey paper of the authors [34].
The result of the convolution is rescaled in amplitude by a non-linear factor:
(3)
where is the signal amplitude obtained from the convolution-autocovariance calculation and is the amplitude of the original signal.-
3.. We apply an envelope detector to the outcome of the IFT (inverse Fourier transform) from the previous step. This demodulates the signal thanks to the loss of symmetry after the amplitude and phase rescaling.
-
4.. Finally, the filtered signal is recovered after a division of the result of the envelope detector.
To summarize, we represent the whole in Figure 2, where the input signal is the acoustic signal obtained directly from the induction motor.
For the sake of completeness, we give more details about using second-order statistics to explain the result of the convolution of a signal with its autocovariance. The autocovariance measures the scattering of the signal around the mean value [33,34,48,49]. If is the signal, then:
(4)
where is the mean and is the autocorrelation of . If is zero then the autocorrelation matches the autocovariance. For ergodic stationary data, we may assume to be constant for all . Then, the autocovariance function of (3) becomes:(5)
2.3. Convolution-Autocovariance Calculation
When acquiring a periodic signal , the noisy observed signal is described as:
(6)
where and are the input signal and an additive stationary noise. Furthermore, , , and for , are the amplitude, frequency, and phase for every single harmony component. Using this representation, we obtain:(7)
The autocorrelation of a white noise with zero mean and variance is a Dirac distribution . If the signal is harmonic, it returns a harmonic signal with null phase. Then from (4) and (6), we obtain:
(8)
From (7), it is noted that the phase information of the original data is missed. For retrieving the phase information we can convolve . The calculation process is described below:
(9)
If the number of data samples tends to infinity, , then we recover the phase information:
(10)
2.4. Spectral Pattern Recognition for Broken Bar Detection
For the sake of clarity, we set a method that enables us to classify damages failures, see Figure 3. It is based on the descending order spectrum applied to the signal obtained after noise reduction. Then, we apply spectral subtraction respect to the healthy motor signal. Finally, we apply a moving average block to smooth the signal by impulsive components of the spectral subtraction. We point out that the pattern recognition algorithm requires the use, as a basic pattern, of the healthy motor signal. This is not a major constraint since this signal can be obtained easily during motor commissioning. At any moment of the posterior motor operation, actual motor samples can be obtained and compared with the healthy one, analogous to what is done in adaptive systems schemes.
3. Results
3.1. Comparison and Assessment of the Proposed Noise Reduction Algorithm: Signal-To-Noise Ratio
We first check the robustness of the proposed algorithm in environments with a variable signal-to-noise ratio. To this end, we compare it with the algorithm proposed in [5]. This works it is based on EMD, which is a non-stationary and non-linear time series analysis method. The method decomposes the series into vectors or intrinsic functions, which are obtained, by decomposing the main function [50,51].
Several decomposition modes were used to obtain the desired spectral components, as shown in [5], see Figure 1 and Table 1. The component at 60 Hz was obtained there with the fifth decomposition mode IMF 5. In the experiments carried out here, we use similar signals (a signal based on four harmonics with frequencies of 60 Hz, 200 Hz, 780 Hz, and 800 Hz, respectively) and similar values of signal-to-noise ratio (−10.9172 dB). The initial correlation value with respect to the harmonic signal without noise was 0.1756, and the final correlation value after processing using the proposed noise reduction algorithm was of 0.7717.
A signal decomposition can bring disadvantages with respect to a method based only on the extraction of statistical characteristics from the signal, since a priori we do not know how many levels of decomposition a deterministic harmonic signal immersed in a random signal can have and neither do we know if each harmonic will have the same decomposition levels.
In the case of the method proposed in [23], various levels of decomposition were used for each harmonic contained in the useful signal. However, the algorithm proposed in our work processes the harmonic signal in its entirety independent of its noise content. Figure 4 and Figure 5 show the signal without noise and the spectrum of the signal contaminated with noise.
In turn, Figure 6 shows a comparison of the result once the harmonic signal taken as a sample has been processed. This permits us to test the influence of the signal–to-noise ratio in the proposed algorithm.
In Figure 6, it can be observed that after processing, the desired spectral components are visible and the noise was reduced. Moreover, with our proposed algorithm, the use of a signal decomposition mode is not necessary, unlike the algorithm proposed in [23]. It only needs the amplitude, frequency, and phase of the input signal.
As can be seen in Figure 6, the amplitude of the signal at the output is affected during the noise-reduction process, and even more if it has a low SNR ratio at the input. However, the amplitude adjustment process can be improved by applying the result shown in Equation (3) as a non-linear amplitude adjustment factor.
Since the statistical characteristics are based on the expected value operator (continuous model), and in practice a finite value of samples is used, the results will be more similar in amplitude if the number of samples to be processed increases (see Figure 7). In addition, the SNR ratio at the input also influences since the method significantly reduces noise, but at the output there are always samples of residual noise that is not eliminated and influence the amplitude of the signal.
In relation to other classic noise reduction methods for the same experimental case and using the same harmonic signal, a wavelet-based method was used. This method uses soft heuristic Stein’s unbiased risk estimate (SURE) thresholding to obtain the actual threshold, to reduce noise. Table 2 summarizes the results obtained in comparing the proposed method with that based on wavelets.
As can be seen in Table 2, one of the most significant disadvantages of wavelet-based methods is a priori ignorance of the calculation of the noise threshold to be eliminated. This calculation is based on knowing the noise power to reduce the random signal, but in real situations, this parameter is generally unknown An example of application is provided in Figure 8.
In general, the choice of one method or another to reduce noise depends mainly on the experimental situation required to process the information to reduce the random signal. There are different methods, each of which has specific characteristics. The effectiveness of one or the other depends on the characteristics of the process. A combination of methods based on statistical characteristics, as in our case, has the advantage of basing their analysis on the analytical foundations of the random process itself.
3.2. Computational Cost
Likewise, the computational cost is related to the number of samples to be processed, that is, the length of the data window, as well as the number of basic operations, quantified in multiplications and accumulations. Figure 9 shows an estimate of the execution time curve of the entire proposed algorithm, in general, including the noise-reduction process and the fault pattern recognition stage. The execution time increases exponentially as the number of samples to be processed increases. This is mainly due to the multiplication and accumulation processes found in the covariance and convolution operations markedly in the noise-reduction and power spectrum calculation processes.
3.3. Failure Detection
To distinguish among the different types of failure, we have taken acoustic samples of a 1.1 kW motor working at full load. The electric motor was a 4-pole machine coupled to a direct current (D.C.) machine acting like a load. In Figure 10, a test-bench and set of tested rotors are shown.
The processed signals of a rotor with 28 bars present the following characteristics: (1) healthy, (2) one broken bar, two broken bars in the relative positions 1 and 2, (3) two broken bars in the relative positions 1 and 3, and (4) two broken bars in the relative positions 1 and 5. See Figure 10 for the T-bench of the tested rotors.
In each test, we recorded the acoustic noise signal with a conventional smartphone equipped with an internal microphone (type omnidirectional condenser) that enabled us to capture the required signals at a sampling rate of 16 kHz. The sampling window consists of 80,000 data measurements. The sampling window consists of 80,000 data measurements. To prevent any vibration influence in the recording, we have located the microphone was at the same place in all the tests. [23]. We removed the undesired components of the spectrum through the noise-reduction algorithm. Figure 11 and Figure 12 show the signal of the healthy rotor before and after the noise reduction with the harmonic peaks.
We show the results of the one broken bar signal. In Figure 13, Figure 14 and Figure 15, we show the signal spectrum before and after the noise reduction. With zoom in the range of 500–1200 Hz we see a third harmonic that has not appeared in the healthy signal in Figure 11. To solidify the finding of the increment of harmonic peaks on the broken bar signals, we have also looked for it in the rest of the two broken bar signals.
We show the spectrum results for the two broken bars signals in Figure 16, Figure 17 and Figure 18. Regarding the spectrum of the signal with two broken bars in the relative position 1–2, we note four harmonic peaks, with different frequencies with respect to the healthy motor spectrum. Here, the most prominent harmonic amplitude is twice the highest of the previous signal (healthy and with one broken bar). The four harmonic pattern is also present in the signal corresponding to the two broken bars signals in relative positions 1–5, see Figure 16. However, in the spectrum of the two broken bars signals in relative positions 1–3, we only appreciate two peaks, see Figure 17.
We note a characteristic component swaying on 750 Hz frequency in all signals. However, we show the differences in the spectrum of the noise-reduced signals. We show a comparison of all signals spectrum after the noise-reduction process in Figure 19. We show the characteristics peaks of the different signals together.
The novelty of the proposed calculation is that it does not require the frequency component to identify the fault. It is based on the stationary or cycle stationary nature of the acoustic signal and only needs to process the amplitude, frequency, and phase parameters. The obtained bar recognition patterns allow us to recognize between the faulty and the healthy motors. They will let us conclude if the rotor is healthy or if it has one or two broken bars, regardless of their relative position, as we will see in the next section.
4. Discussion
To check the findings from the proposed noise-reduction algorithm, we conduct a correlation analysis between the signal without noise and the signal after being processed to reduce the noise. This also lets us know that we only have noise and interferences that do not provide information.
Figure 20 shows the correlation values’ behavior to the output of the proposed noise reduction algorithm as a function of the signal-to-noise ratio (SNR) input values. We see that the correlation values decrease exponentially as the SNR ratio at the input decreases.
With a threshold correlation value of 0.5, we have a linear intrinsic relationship between the input and the output. The algorithm emits correlation values above 0.5 for an input SNR of −15dB (low or ambient with high noise concentration), which is relevant since it does not require additional information apart from the frequency, amplitude, and phase components in the signal to be processed.
In Figure 21 and Figure 22, respectively, we show the spectral pattern and another graphic making zoom, that identifies the number of broken bars shown within the motor. We point out that there is a unique pattern when the motor presents one broken bar and another unique pattern when it has two broken bars that do not depend on the relative position of the bars.
The procedure is based on the spectral subtraction between the healthy motor’s power spectrum and the power spectrum of the damaged motor. The spectral patterns obtained from subtraction converge to zero since the proposed algorithm is based only on highlighting the differences in the power spectrum according to the fault treated and with respect to the spectrum of the healthy motor. We also see that the result is not influenced by the number of samples to be processed, which is 80,000 in this case.
To prove that there is a relation between the spectral patterns with two broken bar, we compute the Pearson relationship coefficients of the vectors (red, green, blue) inferred from the output appeared in Figure 18. The results are presented in Table 3. We see that the Pearson correlation coefficients are close to one, which yields that we are in the presence of a common pattern when the engine has two broken bars independently of the relative of the bars.
5. Conclusions
We have shown that there are common characteristics within the amplitudes and frequencies of signals compared to different kinds of motor with broken bars. We have presented how the acoustic signal of an induction motor is processed. In the first analysis, a noise-reduction process is proposed based on convolving the signal with its covariance function. The adjustment amplitude process in the spectral domain is used to retrieve the original signal without noise.
On the other hand, another algorithm for the identification of broken bars independent of the relative position was proposed. The method was based on the spectral subtraction of every processed signal’s descent-sorted power spectrum. With this algorithm, a common spectral pattern was obtained when the motor presents two broken bars and a unique pattern when the motor has one broken bar.
Results are highly satisfactory, given by the significant signal-to-noise ratio achieved in the noise reduction process and the identification spectral pattern obtained from the broken bars.
Author Contributions
Conceptualization, M.E.I.M., J.A.A.-D., P.F.d.C., and J.A.C.; Methodology, M.E.I.M.; Software, M.E.I.M.; Validation, M.E.I.M., J.A.A.-D., P.F.d.C. and J.A.C.; Formal Analysis, M.E.I.M.; Investigation, M.E.I.M., J.A.A.-D., P.F.d.C. and J.A.C.; Resources, M.E.I.M., J.A.A.-D., P.F.d.C. and J.A.C.; Data Curation, M.E.I.M. and J.A.A.-D.; Writing—Original Draft Preparation, M.E.I.M., J.A.A.-D., P.F.d.C. and J.A.C.; Writing—Review and Editing, M.E.I.M., J.A.A.-D., P.F.d.C. and J.A.C.; Visualization, M.E.I.M.; Supervision, J.A.A.-D., J.A.C. All authors have read and agreed to the published version of the manuscript.
Funding
This research was funded by MEC, grant number MTM 2016-7963-P; Spanish ‘Ministerio de Ciencia Innovación y Universidades’ and FEDER program in the framework of the ‘Proyectos de I+D de Generación de Conocimiento del Programa Estatal de Generación de Conocimiento y Fortalecimiento Científico y Tecnológico del Sistema de I+D+i, Subprograma Estatal de Generación de Conocimiento’ (ref: PGC2018-095747-B-I00); and Generalitat Valenciana, Conselleria de Innovacion, Universidades, Ciencia y Sociedad Digital, (project AICO/019/224).
Conflicts of Interest
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.
Figures and Tables
Figure 1. Illustration of: An example of industrial motor (2000 H.P) with broken rotor bars.
Figure 3. Block diagram of the spectral pattern recognition algorithm for identifying broken bars.
Figure 4. Harmonic signal without noise used to compare our method with that described in [5].
Figure 5. Spectrum of the noisy signal used to compare our method with that described in [5].
Figure 6. The spectrum of the original signal (blue), the spectrum of the noisy signal (black), and the spectrum obtained after applying our algorithm for noise reduction without additional information of the input signal (red); 16,000 processed samples.
Figure 7. The spectrum of the original signal (blue), the spectrum of the noisy signal (black), and the spectrum obtained after applying our algorithm for noise reduction without additional information of the input signal (red); 64,000 processed samples.
Figure 8. The spectrum of the original signal (blue), the spectrum of the noisy signal (black), and the spectrum obtained after applying the wavelet algorithm for noise reduction.
Figure 9. Execution time in terms of the number of samples. The execution time was calculated using the Matlab tictoc function and on a computer with an i5 processor and 4Gb of RAM.
Figure 11. Healthy rotor spectrum before the noise reduction with the harmonic peaks.
Figure 12. Healthy rotor spectrum after the noise reduction with the harmonic peaks.
Figure 13. One broken bar signal spectrum before the noise reduction with the harmonic peaks.
Figure 14. One broken bar signal spectrum after the noise reduction with the harmonic peaks.
Figure 15. One broken bar signal spectrum before the noise reduction with zoom in the range of 500–1200 Hz.
Figure 16. Spectrum of the two broken bars, in relative position 1–2, signal after the noise reduction with the harmonic peaks.
Figure 17. Spectrum of the two broken bars, in relative position 1–5, signal after the noise reduction with the harmonic peaks.
Figure 18. Spectrum of the two broken bars, in relative position 1–3, signal after the noise reduction with the harmonic peaks.
Figure 19. A comparison for all signals spectra after the noise-reduction process.
Figure 20. Output correlations values in terms of the input signal to noise ratio (SNR).
Figure 21. Spectral pattern obtained from the application of the proposed pattern recognition algorithm showed in Figure 3, to identify broken bars.
Figure 22. Zoom of the spectral pattern obtained from the application of the proposed pattern recognition algorithm showed in Figure 3, to identify broken bars.
Comparison of our proposed algorithm with the competing method proposed in [23].
Proposed Algorithm | Competing Method Ref. [23] | |
---|---|---|
System characteristic | One Input/One Output | One Input/One Output |
Signal to Noise Ratio Input | −10.9172 | −10.9172 |
Input Signal | Harmonic Signal with component (60 Hz, 200 Hz, 780 Hz, and 800 Hz,) | Harmonic Signal with component (60 Hz, 200 Hz, 780 Hz, and 800 Hz,) |
Signal Decomposition | NO | Yes |
Summary of the characteristic for the comparison with the method using wavelet function.
Proposed Algorithm | Competing Method Based on Wavelet | |
---|---|---|
System characteristic | One Input/One Output | One Input/One Output |
Signal to Noise Ratio Input | −10.9172 | −10.9172 |
Input Signal | Harmonic Signal with component (60 Hz, 200 Hz, 780 Hz, and 800 Hz,) | Harmonic Signal with component (60 Hz, 200 Hz, 780 Hz, and 800 Hz,) |
Signal Decomposition | NO | NO |
Thresholding Calculation | NO | Yes |
Obtained Correlation | 0.7717. | 0.3303 |
Pearson correlation coefficient values between output signals with two broken bars.
Signal | Correlation Value |
---|---|
Bar 1–2 vs. Bar 1–3 | 0.9883 |
Bar 1–2 vs. Bar 1–5 | 0.9690 |
Bar 1–2 vs. Bar 1–5 | 0.9812 |
You have requested "on-the-fly" machine translation of selected content from our databases. This functionality is provided solely for your convenience and is in no way intended to replace human translation. Show full disclaimer
Neither ProQuest nor its licensors make any representations or warranties with respect to the translations. The translations are automatically generated "AS IS" and "AS AVAILABLE" and are not retained in our systems. PROQUEST AND ITS LICENSORS SPECIFICALLY DISCLAIM ANY AND ALL EXPRESS OR IMPLIED WARRANTIES, INCLUDING WITHOUT LIMITATION, ANY WARRANTIES FOR AVAILABILITY, ACCURACY, TIMELINESS, COMPLETENESS, NON-INFRINGMENT, MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. Your use of the translations is subject to all use restrictions contained in your Electronic Products License Agreement and by using the translation functionality you agree to forgo any and all claims against ProQuest or its licensors for your use of the translation functionality and any output derived there from. Hide full disclaimer
© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.
Abstract
Featured Application
We provide a non-intrusive tool for the detection of adjacent and non-adjacent bar breakage from the acoustic noise radiated by a motor. It can be included as a smart application in a transportable device.
AbstractWe apply power spectral analysis based on covariance function and spectral subtraction to detect adjacent and non-adjacent bar breakages. We obtain a spectral pattern when the signal presents one or various broken bars, independent of the relative position of the bar breakages. The proposed algorithm gives satisfactory results for detectability compared to some previous research. Additionally, we also present illustrations of faults and signal to noise in the noise-reduction stage.
You have requested "on-the-fly" machine translation of selected content from our databases. This functionality is provided solely for your convenience and is in no way intended to replace human translation. Show full disclaimer
Neither ProQuest nor its licensors make any representations or warranties with respect to the translations. The translations are automatically generated "AS IS" and "AS AVAILABLE" and are not retained in our systems. PROQUEST AND ITS LICENSORS SPECIFICALLY DISCLAIM ANY AND ALL EXPRESS OR IMPLIED WARRANTIES, INCLUDING WITHOUT LIMITATION, ANY WARRANTIES FOR AVAILABILITY, ACCURACY, TIMELINESS, COMPLETENESS, NON-INFRINGMENT, MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. Your use of the translations is subject to all use restrictions contained in your Electronic Products License Agreement and by using the translation functionality you agree to forgo any and all claims against ProQuest or its licensors for your use of the translation functionality and any output derived there from. Hide full disclaimer
Details




1 Departamento de Telecomunicaciones, Universidad de Pinar del Río, Martí #270, Pinar del Río 20100, Cuba;
2 Instituto Universitario de Matemática Pura y Aplicada, Universitat Politècnica de València (UPV), Camino de Vera s/n, 46022 Valencia, Spain;
3 Instituto Tecnológico de la Energía, Universitat Politècnica de València (UPV), Camino de Vera s/n, 46022 Valencia, Spain