Robust Switching Control and Subspace Identification for Flutter of Flexible Wing
Active flutter suppression and subspace identification for a flexible wing model using micro fiber composite actuator were experimentally studied in a low speed wind tunnel. NACA0006 thin airfoil model was used for the experimental object to verify the performance of identification algorithm and designed controller. The equation of the fluid, vibration, and piezoelectric coupled motion was theoretically analyzed and experimentally identified under the open-loop and closed-loop condition by subspace method for controller design. A robust pole placement algorithm in terms of linear matrix inequality that accommodates the model uncertainty caused by identification deviation and flow speed variation was utilized to stabilize the divergent aeroelastic system. For further enlarging the flutter envelope, additional controllers were designed subject to the models beyond the flutter speed. Wind speed was measured online as the decision parameter of switching between the controllers. To ensure the stability of arbitrary switching, Common Lyapunov function method was applied to design the robust pole placement controllers for different models to ensure that the closed-loop system shared a common Lyapunov function. Wind tunnel result showed that the designed controllers could stabilize the time varying aeroelastic system over a wide range under arbitrary switching.
Aeroelastic problem has been always attracting extensive attentions in past few decades. Along with the rapid development of aeronautical science, the requirements of endurance and maneuverability are constantly growing, leading to lightweight and flexible characteristics so that the aeroelastic problem is more noticeable. Flutter is a self-excited phenomenon that occurs when fluid couples energy into structure motion causing instability and excessive vibration. So flutter suppression is one of the major objective problems in aeroelastic control.
Over the past decades, many research results have been obtained by using different actuator and algorithm to suppress flutter. In early 1980s, NASA Landley center launched active flexible wing plan (AFW) and Benchmark technology projects to suppress the wing model flutter under subsonic and transonic case by using flap [1, 2]. In recent years, a series of studies are conducted to show the efficiency of flap to suppress the bending torsion flutter and limit cycle oscillation [3–5].
Piezoelectric material is also a major topic of discussion due to its advance of broadband frequency range and ease to be embedded. Guan et al.  conduct a flutter suppression experiment to the lifting surface by using distributed piezoelectric actuator. Since the driving force of piezoelectric ceramic is relatively small, advanced actuators derived from piezoelectric material are applied to the aeroelastic control. Ardelean et al.  designed a V-stack piezoelectric actuator to suppress the subsonic flutter of wing section which is capable of producing large force and stroke to drive the flap over a wide frequency range. Burnham and Pitt  conduct the research on buffeting alleviation for the vertical tail by using the macro fiber composite actuator, which is proved to be sufficiently forceful and less influential to the structure.
In addition to the experimental research, control algorithm is also playing an important role in active flutter suppression. Due to the variation of flying condition and model error, uncertainty of system model is inevitable. In order to cope with the impact of uncertain, robust control algorithms are applied to ensure the stability and performance of the closed-loop system. Han et al. thought over the parameter uncertainty and structure uncertainty of aeroelastic system and designed the LQG and controller synthesis, respectively, by comparing two schemes implemented on a vertical thin plate, exhibiting a better performance and scope of action to the flutter control which improve the closed-loop critical speed . Zeng et al. designed a controller for a semifuselage wing subjecting to the model uncertain and gust distribution based on ARMA model and the dynamic pressure could be improved in the presence of uncertainty . Blue and Balas construct the LPV model by taking the dynamic pressure and Mach number as parameter perturbation, by separating the nominal model and uncertain model using linear fractional transform; the LPV controller was designed to suppress flutter of time varying system . A unique controller could only be effective in a limit range. To enlarge the flutter envelope, Huang et al. employed the closed-loop identification to obtain the dynamic models beyond the critical speed and robust controllers are synthesized based on the corresponding model to suppress the wing model within the range between two identified models . McEver et al. conducted parameterization method which make the frequency range selectable to identify the closed-loop model and suppressed the flutter by PPF controller .
Although the corresponding controllers designed for different models are mentioned above, the switching process between each controller has not been fully considered. Even if every controller ensures the performance under the static condition, the stability and performance under arbitrary switching could not been guaranteed . The system stability is ensured by satisfying the inequality according to the theory of Lyapunov stability of linear system, that educes a Lyapunov function to be negative definite. However, for the switching system, the dynamic behavior is substantially nonlinear, which cannot be analyzed by the primitive theory of Lyapunov stability. The Lyapunov function of the switching system has a risk to increase at the time of switching. Therefore, it is essential to study the active flutter suppression under switching process.
The stability of arbitrary switching control is so far attributed to the common Lyapunov function method [15, 16], multiple Lyapunov function method [17, 18], and dwell time method [19, 20]. The common Lyapunov method ensures that if each subsystem shares a common Lyapunov function, then the switching system could maintain the stability under arbitrary switching. Switching synthesis considers the design of controllers as a unified process, each controller is designed relatively to ensure the performance under random switching. The common Lyapunov function theory has been used for analyzing the individual linear model and partial linear model such as fuzzy T-S system . Since the common Lyapunov function based controller design is essentially conservative , multiple Lyapunov function method is proposed. The switching stability can be ensured even if each subsystem has an individual Lyapunov function under certain conditions . Dwell time method is emphasis on the critical switching time for system stability. It can be precisely computed based on exact dynamic model and the conservative of controller design is significantly decreased while the switching speed is less than dwell time. For slow switching system, several efficient methods are proposed based on dwell time theory [22–24].
In this study, system analysis and controller design are all based on the common Lyapunov theory since the dynamic model of the aeroelastic system is time varying and uncertain so that the certain conditions based on exact models are not available. The contents of this paper are as follows: dynamic model of flutter is obtained by system identification, and then multicontrollers are designed by robust pole placement. Dynamic model beyond the critical wind speed could be obtained by using frequency domain subspace method under closed-loop conditions. The common Lyapunov function and liner matrix inequality are used to correlate with each controller to ensure the arbitrary switching stability. Wind tunnel test is implemented to indicate that the flutter of the wing under a specific wind speed can be suppressed by each robust controller and stabilization of arbitrary switching among controllers can be ensured perfectly. The main contributions of this research are that the switching stability is completely insured subjecting to the wind speed variation and the robust switching control method is effective significantly in suppressing the wing flutter.
2. Wind Tunnel Model and Aeroelastic Equation
2.1. Model Description
A wing model is designed to verify the control algorithm in wind tunnel test. The model is made of elastoplastic with elastic modulus 1.1 Gpa, density 1430 , and Poisson ratio 0.39 and the cross section of the wing model is NACA0006. The max thickness is 6% of the chord length. A lumped mass of 36 grams is mounted on the tip of the wing in order to reduce the flutter speed in wind tunnel test. The wing root is vertically fixed on the base of the closed wind tunnel. Real model and low speed wind tunnel are shown in Figure 1. Due to the limitation of maximum wind speed, wind tunnel test of flutter suppression is carried out within the flow speed range below 35 m/s.
The configurations of actuator and sensor are shown in Figure 2. MFC is adhered to the root of the wing model as the control actuator and accelerometer are mounted on the tip of the wing as the monitoring and feedback signal.
Modal parameter of the wing can be obtained via FEM analysis or system identification test. The modal analysis is implemented in Nastran FEM software subjected to the first four modes. In system identification test, a low-pass white noise is applied to the MFC as input and the accelerometer signal is measured as output. The natural frequency of the wing model can be identified from the input and output signals via subspace method that will be mentioned later. Comparison between the numerical and experimental results is shown in Table 1.
The difference between the numerical and experimental results is contributed to the deviation of the material parameters and the neglection of MFC actuation. Analysis and test results later show that the flutter occurs in the third mode at wind speed 28 m/s.
2.2. Dynamic Equation of Aeroelastic System
MFC actuator shown in Figure 3 is internally arranged by a combination of piezoelectric fibers, which can be applied to the piezoelectric material by means of a cross electrode to produce a large driving force. The constitutive equations of piezoelectric materials can be expressed as follows:where , , , are the strain, stress, electric displacement, and electric field, respectively, and , , denote elastic compliance and piezoelectric constant.
The equivalent stress of the driven force as a result of applied voltage on the cross electrode can be expressed as follows:And the moment applied by the MFC to the wing model can be expressed as follows:where denote the distance between stress point and middle layer of the airfoil. The integral domain includes all the areas covered by the piezoelectric material.
To determine the aeroelastic equation, Hamilton principle is used derived from the virtual work of potential energy, kinetic energy, and generalized external force. The total potential energy and kinetic energy of the aeroelastic model can be given bywhere is density variable and is displacement variable. The generalized external force of the system is composed of MFC actuator force and aerodynamic force, which can be expressed aswhere is unsteady aerodynamic matrix which can be computed via doublet lattice method . are reducing frequency and Mach number, respectively. By applying Hamilton principle, the variation formulation can be established by
In order to obtain the equations of motion in modal space, coordinate transformation of the displacement variable is carried out. Displacement of the wing model is expressed in terms of generalized coordinates:Substituting (7) into (6) and performing the variation operation in terms of , expression of the aeroelastic equation can be obtained aswhere , , are modal mass, modal damping, and modal stiffness matrix, namely,and is the unsteady aerodynamic matrix under modal coordinate and is the dynamic system from the control output to the generalized force applied to modal variable.
2.3. Identification of Aeroelastic System
In this study, frequency domain subspace identification method is used to obtain the open-loop and closed-loop aeroelastic system. The identification algorithm is implemented based on the frequency response data to estimate system matrix of state space model and it exhibits a reliable performance for open-loop system and closed-loop system identification .
Computation simplicity is the major advantage over other method. Frequency domain estimation ensures that the identified data be selectable near the flutter frequency that provides an accurate information. Sufficient accuracy can be obtained if the outside disturbance is small enough. The state space description of multi-input and multi-output system in frequency domain can be expressed as:where , , are the discrete Fourier transformations of , , . In this study, the aim of system identification is to obtain the system matrices , , , . is the response of signal by the acceleration transducer adhered on the wing model. is the output signal of the controller. is the generalized operator. For the frequency point in frequency domain, denotes the corresponding value of circular frequency. is sampling frequency. By iteration of (10), the high order equation can be obtained byRewrite into matrix form:where is the generalized observable matrix, which can be expressed as is the lower triangle Toeplitz matrix:Several values in the frequency domain are selected and all satisfied (12), rewriting into matrix form:where , is the number of discrete points in frequency domain, and , are the Van Dermond block matrix expressed as:Separating (15) into real part and imaginary part, it can be expressed aswhere ,
To estimate , decomposition is implemented to the matrix of input and output matrix in frequency domain:The singular value decomposition is implemented to matrix :By applying orthogonal projection to the row space of , Consequently we can get the estimation values of asAccording to (21),The least square method was utilized to estimate and ; thus, namely, The transfer function of the control system from to isand it can be written in the form of Kronecker productwhere represents a column vector of the matrix which is arranged in order and then and can be obtained by using least square method.
3. Controller Design
Flutter is caused by the instability of the aeroelastic system in general and the poles of the system should be located in the left half of complex plane for the linear stability system. Now, pole placement controllers are designed to make the closed-loop system poles in the left half of the complex plane beyond the flutter speed. Since aeroelastic system is time varying according to the change of the Mach number, attack angle, dynamic pressure, and other factors, single controllers cannot match well with the model and would be invalid if the factors vary in large scope. So single controller is limited in suppressing the wing flutter. In this research, several robust pole placement controllers are designed based on the aeroelastic model under several specific wind speed range; therefore each controller can theoretically place the closed-loop system poles within the designated area in complex plane.
After the state space model of aeroelastic system is established, linear fractional model of uncertain system shown in Figure 4 is taken into consideration.
3.1. System Modeling
The linear fractional model in form of additive uncertainty is consider as a feed forward dynamic uncertainty which is added to the model that represents the dynamic deviation caused by the model identification and wind speed variation. The state space with additive uncertainty can be expressed aswhere represents an uncertain complex gain matrix satisfying and is a prescribed constant which can be explained as the peak error of the model in frequency domain. , , , are the state space realization of , respectively, shown in Figure 4. As the uncertainty is represented by the form of additive, it can be inferred that .
Consider a dynamic feedback controller which can be expressed asSubstituting (27) into (26), the closed-loop system can be expressed aswhere the closed-loop system of uncertain model can be expressed asThe purpose of the controller design is to place the poles of the uncertain closed-loop system within a specific region.
3.2. Robust Pole Placement Controller Design
LMI Region. A LMI region is any subset of the complex plane with characteristic function:where , are the Hermitian matrices. For the region of left half plane and , the characteristic function can be expressed asThe symbol “<0” denotes that is negative definite. Following content introduces the detail of realizing (31) for the closed-loop system.
Definition 1. The uncertain system (27) is robustly -stable if the eigenvalues of lie in for all admissible uncertainties .
Lemma 2 (, robust -stability). The sufficient and necessary conditions for the robust -stability of the closed-loop system are that a positive definite symmetric matrix which satisfies with (32) is existence:Equation (32) is a nonlinear matrix inequality for ; by using variable substitution, it can be converted to the linear form.
Theorem 3. A output feedback controller and a symmetric matrix are available such that (32) holds if and only if two symmetric matrices and state space matrices are existence, namely,
Proof. The symmetric positive definite matrix that satisfied (32) can be expressed as the following form:Two matrices and are defined byAnd it can be verified thatIt has been proven  that the full rank matrices , are necessarily existence if satisfied (32), then it can be inferred that and are invertible. Pre- and postmultiplied by and , this yields Define the variable substitution equation:Pre-and postmultiply (32) by and and pre-and postmultiply by and , substituting (39) and into the result, then (33) and (34) are satisfied:Since , , and the necessity is proved.
Similarly, pre- and postmultiply (33) by and and pre- and postmultiply (34) by and : If the above LMIs are feasible, the matrices , , are derived as follows:1.Calculate the matrices , based on .2.Solve (39) for , , , .
3.3. Switching Controller Design
Equations (33), (34), and (39) give the synthesis process of robust controller. However, the robust controller will lose efficiency if the system uncertainty is larger than . A group of switching stabilized controllers are synthesized according to the time varying aeroelastic model and its additive uncertainty under different wind speed, shown in Figure 5.
The uncertain switched system of aeroelastic model in state space form can be expressed as follows:where is the time varying switching signal in set that represents different model under the time variant condition. In this study, wind speed is considered as the characteristic condition which arbitrarily varied with time. The switching signal is determined by decision variable which reflect the varying of wind speed. Each subsystem of (43) represents the state space model of different wind speed. A set of controllers are designed by the following theory to ensure the stability during arbitrary process.
Theorem 4. System equations (42) have a common Lyapunov function if there exists .
Proof. It can be seen that if the closed-loop system (42) is satisfied with (43) and (44), there exists a common and which satisfy with (43). For each subsystem, pre- and postmultiply by and , namely,From Lemma 2, (47) is equivalent to the inequality which indicated that each subsystem shares a common Lyapunov function.
4. Wind Tunnel Test of the Wing Flutter Suppression
4.1. Description of Experimental Setup of Wind Tunnel Test
The aeroelastic control system, shown in Figure 6, consists of embedded controller, accelerometer, wind speed transducer, signal conditioner, power amplifier, and MFC actuator. The acceleration response signal is acquired by the sensor and fed back to the controller. NI PXIe-1071 controller is used as the measured signal processing and control algorithm implementation. PXIe-6356 16 bit input and output module which is integrated with the controller performs a high speed and precision measurement that ensure supporting the control performance. Real-time performance of the control system can be obtained by the NI LabVIEW software and high calculation of the controller. The control and signal processing algorithm are programmed on LabVIEW language. Control signal derived from the control algorithm operation is applied to MFC actuator via the power amplifier. The controller is designed at wind speeds 27, 31, and 33 m/s, respectively, and discretized at the sampling time 0.002 s. Experimental data are recorded and displayed on the platform to verify the control performance.
In order to verify the dynamic characteristic and stability under arbitrary switching of the closed-loop system, wind speed is manually regulated by the knob on the control cabinet of wind tunnel, shown in Figure 7.
(a) Experimental apparatus of wind tunnel test
(b) Control cabinet of wind tunnel speed
4.2. Model Identification of Aeroelastic System in Wind Tunnel
From Section 2.3 we know that the state space model of aeroelastic system can be obtained by frequency domain subspace method. A 25 Hz to 45 Hz band pass signal was used as the input signal which is applied to the MFC actuator via power amplifier and the accelerometer response signal is used as the output signal of the model. In this study, frequency response function of the input and output signals is estimated as the identified data by the ratio of cross power spectrum and self-power spectrum as follows:
Power spectrum can be estimated based on the time domain data via Welch method. Three frequency response models are estimated under the condition with wind speeds 27, 31, and 33 m/s, respectively. The identification results of three frequency response function of input and output signals are shown in Figure 8.
(a) FRF under the condition with 27 m/s wind speed
(b) FRF under the condition with 31 m/s wind speed
(c) FRF under the condition with 33 m/s wind speed
From Figure 8, we find that the major energy of flutter is concentrated near the range of flutter frequency. To simplify the process of the controller model, all modes except for the flutter mode will be neglected in identification process. A second-order state space model is used to approximate the aeroelastic system by mean of optimal fitting in frequency domain. The identified models in transfer function forms of three different wind speeds by subspace method are as follows:
From the comparison of frequency response function between identified and test models, we find that the modulus error and phase error caused by identification deviation and variation of flow speed are existence. It will be compensated in the process of robust controller design.
4.3. Robust Pole Placement Control Test in Wind Tunnel
Vibration divergence of the wing model occurs when wind speed exceeds the critical value shown in Figure 9; however, it would not increase infinitely. Due to hardening characteristics of the structure, approximate stable vibration appears after a specific magnitude is attained.
The robust pole placement controllers are designed offline based on above three identified models, and each controller is solved by LMI.
By assigning the pole position to a specific region and estimating the uncertainty according to the frequency response, the robust pole placement controller under the condition witch 27 m/s, 31 m/s and 33 m/s is designed by setting
Three identified state space models with different wind speed are as follows:
Before flutter suppression experiment, the designed controller should be discretized. Zero-order hold transformation is used to convert the analog controller to the form of difference equation. The discrete transfer functions of designed controller in three wind speed are
To demonstrate the effectiveness of designed controller, three tests are implemented with 28 m/s, 31 m/s, and 33 m/s single wind speed. During the test, controllers are switched off at first and then vibration is converged quickly after the controllers were switched on, shown in Figure 10.
(a) Time domain response with 28 m/s wind speed
(b) Time domain response with 31 m/s wind speed
(c) Time domain response with 33 m/s wind speed
The closed-loop frequency response functions are shown in Figure 11.
(a) FRF under the condition with 28 m/s wind speed
(b) FRF under the condition with 31 m/s wind speed
(c) FRF under the condition with 33 m/s wind speed
The comparisons between the power spectrum of the controlled and uncontrolled response signal illustrate that the major vibration energy of flutter and turbulence is efficiently alleviated by the robust controllers shown in Figure 12.
(a) Power spectrums with 28 m/s wind speed
(b) Power spectrums with 28 m/s wind speed
(c) Power spectrums with 33 m/s wind speed
Figure 13 shows the root locus of open-loop and closed-loop aeroelastic system with the variation of wind speed identified by subspace method. We find that the aeroelastic system is stabilized by the controller over the flutter speed and all poles of closed-loop system located in stabilization area.
4.4. Switching Control Test in Wind Tunnel
The switching control signal is obtained by the wind pressure transducer placed on the wall of closed tunnel. Switching rule is expressed by the following expression:
By setting the parameters , , , , , of the switching control and solving the LMI of each subsystem simultaneously, three controllers can share a common Lyapunov function. Then, discretized transfer functions of switching controllers are as follows:By comparison of the damping ratios between open loop and closed loop, it shows that the damping ratio of aeroelastic system is improved by the switching controllers.
The wind speed of wind tunnel is varied from the range of 28 m/s to 35 m/s which was controlled by the knob of the panel. The controller is switched dynamically over three candidate controllers and the time domain of the switching signal is shown in Figure 15.
Since the wind pressure signal was disturbed by high frequency noise generated by wind tunnel, there is a serious risk that the switching signal could be changed quickly. The existence of common Lyapunov function for each subsystem ensures the unconditional stability of arbitrary switching. From Figure 16 we find that the switching signal keeps stability while the wind speed changes from 27 m/s to 35 m/s. The vibration amplitude is suppressed about 80%.
Subspace identification is used to model the aeroelastic system in this research. Aeroelastic model under different wind speeds 27 m/s, 31 m/s, and 33 m/s is identified though experiments. The model under 27 m/s is below the critical speed and controller designed based on this model can efficiently enlarge the critical speed and increase the Signal-to-Noise ratio for the identification under 31 m/s, and it has the same effect under 31 m/s. When the flow speed is closed to the critical speed, flutter mode plays a dominant role of vibration, so a two-order model is used to approximate the unstable mode. The instability of vibration can be represented by “negative damping” in vibration. The characteristic equation of second-order transform function is expressed as ; as the damping coefficient is negative, the system is unstable. The identification results in Figure 14 show that the root locus is always in the left half part of the complex plane below the critical speed, the damping coefficient is changed from 0.13 to 0.01 and will always be positive. While the flow speed is higher than the critical speed, the root locus locates on the right half part and the damping coefficient is negative. The identification result shows an accurate reflection to the aeroelastic instability.
The uncertainty of the system is caused by two aspects: change of wind speed and the error of identification. The norm of the uncertain dynamic for controller design is based on the measured frequency response function. The modulus deviation of the FRF under 27 and 31 is less than 0.6, so the norm can be considered as less than this value while wind speed varies from 27 to 31 m/s. Similarly, norm of uncertainty under the wind speed variation from 31 to 33 m/s is set to be 0.8.
In practice, if each controller is individually designed, a global Lyapunov function would not be shared by each closed-loop system so that the convergence can only be satisfied under no switching or slow switching condition. At switching time, total energy of the system is still possible to increase. If the flow speed variation and switching time satisfy a specific condition, vibration is still possible to divergence. The major advantage of common Lyapunov function method is absolutely alleviating the risk of switching instability and ensures the safety of flutter control.
In this paper, macro fiber composite is used to suppress the flutter of wing model under the variation of wind speed. The frequency subspace identification and robust pole placement switching controller design method are discussed and the efficiency of this method is validated though flutter suppression test. The main contributions of the paper are the following two aspects:(1)The subspace identification method in frequency domain is utilized to obtain the dynamic model under open and closed-loop condition. Identified model effectively reflects the instability of aeroelastic system. The robust pole placement controllers are designed based on the identified model and solved by the LMIs. Theoretical analysis and wind tunnel test are implemented to prove that the time varying system within a range of uncertainty can be stabilized by the designed controllers.(2)To enlarge the stability envelope, multiple controllers are applied to suppress the flutter under different wind speed range. Three switching controllers are designed by the common Lyapunov method which ensure that the closed-loop system is stable under arbitrary switching. Theoretical analysis is implemented to prove that the designed controller can stabilize the time varying system within a certain range of uncertainty. To verify the effective of control algorithm, the switching control test is implemented in wind tunnel. Wind speed is manually adjusted to simulate the random change of flow. Experimental result shows that the switching control algorithm can successfully suppress the flutter over a wide range of speed under the arbitrary switching.
|:||System matrix of the state space model|
|:||Input matrix of the state space model|
|:||Output matrix of the state space model|
|:||Feedthrough matrix of the state space model|
|:||State vector of the state space model|
|:||Control signal of aeroelastic system|
|:||coordinate of wing model|
|:||coordinate of wing model|
|:||coordinate of wing model|
|:||Displacement of the wing model|
|:||Symbol of variation|
|:||Output vector of the state space model|
|:||Abbreviation of symmetric matrix|
|:||Max singular value of system matrix.|
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
This work is supported by funding of Jiangsu Innovation Program for Graduate Education, Grant no. KYLX15_0237 (the Fundamental Research Funds for the Central Universities) and Priority Academic Program Development of Jiangsu Higher Education Institutions.
J. Sun and M. Li, “Effects of motion-induced aerodynamic force on the performance of active buffeting control,” Journal of Vibroengineering, vol. 18, no. 2, pp. 951–964, 2016.View at: Google Scholar
D. Guan, W. M Chen, and M. Li, “Flutter suppression using distributed piezoelectric actuators,” Chinese Journal of Aeronautics, vol. 13, no. 4, pp. 211–215, 2000.View at: Google Scholar
P. A. Blue and G. J. Balas, “Linear parameter varying control for active flutter suppression,” in Proceedings of the Guidance, Navigation, and Control Conference, 1997, pp. 1073–1079, usa, August 1997.View at: Google Scholar
X. Zhao, P. Shi, Y. Yin, and S. . Nguang, “New results on stability of slowly switched systems: a multiple discontinuous Lyapunov function approach,” Institute of Electrical and Electronics Engineers Transactions on Automatic Control, vol. 62, no. 7, pp. 3502–3509, 2017.View at: Publisher Site | Google Scholar | MathSciNet