Abstract

Vibration control of wind turbine blade with minimum control cost is investigated. Realization of minimum control cost is based on pole placement with minimum-order observer (PPMO). The blade analysis employs a novel compromise method between the 2D airfoil analysis method and the 3D coupled blade body analysis method based on data fitting. It not only ensures certain accuracy, but also greatly improves the speed of calculation. The Wilson method, developed on the basis of the blade momentum theory, is adopted to optimize the structural parameters of the blade, with all parameters fitted as general model Sin6 (Sum of Sine) fitting curves. Also the aerodynamic coefficients based on data obtained by Xfoil software are fitted. Pole placement technology based on minimum-order observer is applied to control unstable vibrations of vertical bending and lateral bending with minimum control cost characterized by the energy consumption of the controller. The pole placement technology is a novel pole assignment technique based on self-poles derived from constant stable eigenvalues, which can effectively avoid the mismatch problems caused by pole selection. The superiority of PPMO can be apparently demonstrated by comparison of linear quadratic regulator (LQR). Analytical proof of the control accuracy and feasibility analysis of the physical realization of the PPMO algorithm are also investigated by experimental platform of hardware-in-the-loop simulation.

1. Introduction

The flutter analyses of wind turbine blades usually include two cases: typical blade section analysis and the three-dimensional coupled blade body analysis. The former is widely used in quick study of flap/lag vibrations or bending/twist vibrations of single cross-section with a wide range of structural and external parameters [1, 2]. An aeroservoelastic model was presented to study the microtab control of the blade section undergoing moderate stall flutter and deep stall flutter separately [1]. Effects of aeroservoelastic pitch control for flutter suppression on wind turbine blade section subjected to combined flap/lag motions were studied, which was dedicated to solving destructive flapwise and edgewise instability of stall-induced flutter [2]. Although this method is simple and convenient, it cannot fully reflect the vibration behavior of the whole blade body. The latter method involves the study of dynamic analysis of coupled displacements and aeroelastic characteristics of the whole blade body [36], which can more accurately characterize the dynamics of the blade. Li et al. [3] presented the analysis of dynamic characteristics of horizontal axis wind turbine blade, where the mode coupling among axial extension, flap vibration, lead-lag vibration, and torsion was emphasized. Lee et al. [4] investigated the performance and aeroelastic characteristics of wind turbine blades based on flexible multibody dynamics, a new aerodynamic model, and the fluid-structure interaction approach. Members of our project group [5] also investigated the stall-induced aeroelastic stability of wind turbine blade with bending-bending-twist coupling. Andersen et al. [6] investigated the fatigue load reduction applying trailing-edge flaps on wind turbine blade that is another structural form of the blade, rather than the general blade structure. Load control of wind turbine was analyzed, with requirements for load monitoring sensors, sensing technologies, and applications discussed in [7]. However, all this literature involves complicated blade structures and aeroelastic behavior, in which their solutions involve both complex decoupling problems and decomposition of aerodynamic forces along the blade spanwise direction. Whether it is based on numerical analysis or by finite element simulation, the solution process must be a time-consuming process [5]. Needless to say, actual vibration control on the basis of the actual controller hardware is impossible due to the complexity of these algorithms or complicacy of the control system hardware itself.

Structural failures observed in recent years of large wind turbine blades have been mostly attributed to severe vibrations in flapwise and edgewise directions [8]. Some earlier studies involved the flap/lead-lag stability analysis based on eigenvalue problem [9] and stability control based on fuzzy PID algorithm and LQR controller [10], with the actual analyzed structure being single blade section. Members of our group investigated a variety of methods of intelligent control for flap/lead-lag flutter, also based on single blade section [2]. As for analysis of the whole blade body as shown in [3], the numerical difference methods adopted to solve the system integrodifferential equations of the flap vibration and lead-lag vibration are so complex that it is impossible to realize the relevant control algorithms in the actual hardware.

Since the cross-sections of different sizes experience different loads, the combined effects of different loads have far surpassed those of the same loads of the typical sections of the same sizes described in related references [1, 2, 9, 10]. Therefore, the effectiveness of design parameters of cross-sections of different sizes needs to be confirmed by aeroelastic stability analysis for safety reasons. At the same time, the control cost of flap/lead-lag vibrations and the hardware implementation of the control algorithm are not mentioned in above related researches. The “control cost” here mainly refers to the amplitude of the controller, the energy consumed by the controller, or the power required to drive the controller.

In present work, a novel compromise method (NCM) between the 2D typical section analysis method and the 3D coupled blade body analysis method is investigated to study vertical and lateral vibrations of blade tip. Then the analysis accuracy might be between the accuracies of the two structural model analysis methods due to the approximate compromise behavior. The NCM not only preserves the basic characteristics of typical section model, but also extends the integral behavior of blade body along the spanwise direction but avoids the decoupling problem of partial differential equations. More importantly, the vibration control of blade model based on NCM is more efficient and easier to be implemented in controller hardware. In addition, realization of minimum control cost is based on pole placement with minimum-order observer (PPMO), which can be depicted as pole placement technology based on minimum-order observer. The PPMO is applied to control unstable vibrations of vertical/lateral directions based on self-poles derived from constant stable eigenvalue analysis, which is exactly the advantage of the method NCM. At the same time, the minimum-order observer is applied in PPMO to estimate the unmeasurable variables.

It should be stated that the emphasis of present study is on the application of the control method, but not on the aerodynamic nonlinearities associated with dynamic stall in wind turbines [1]. However, the control algorithm proposed here can also provide a method to solve the problem after linearization of the stall nonlinear system.

2. Theoretical Model

2.1. Equations of Motions of Blade Section

Consider the blade structure and coordinate system in Figure 1, the length of the blade is  m in direction and it has the rotating speed , normal to the plane of rotation. The origin of the rotating axis system is located at the rigid root. The characteristic cross-sectional dimension, chord length is at a distance from the hub . It is assumed that the blade section located at the distance is connected by in lead-lag direction and by in flapwise direction , where both and are rigidly connected to the hub . The wind velocity is denoted by . The attack angle is denoted by , with pretwist angle by . It is assumed that the rotating plane will not change due to yaw, and as the design of the pretwist angle increases the torsional stiffness, the elastic torsional displacements of different blade sections can be approximately neglected. The rotating speed can be determined by , according to subsequent structural parameters in Section 2.2.

Consider the single blade section in Figure 1; the equations of motions in and directions are derived using Lagrange’s equations as depicted in [2]. Meanwhile considering dynamic stiffening effect in [9], the equations of vertical () and lateral () motions of blade section are expressed as follows:where the incoming velocity is , with incoming angle being ; the mass per length of blade section is described by ; the natural angular frequencies of motions in and directions are described by and , respectively; the damping ratios in and directions are denoted by and , respectively; the Lift coefficient and Drag coefficient are described as and , respectively.

2.2. Optimal Design of Airfoil Parameters

As the actual sectional position is different, the structural parameters of airfoil and blade body are also changed. In order to take full account of the stability performance, the blade structure is optimized by Wilson method [11, 12], and the aeroelastic model and the aerodynamic coefficients are fitted based on the blade element momentum (BEM) theory [13, 14]. The design uses 21 cross-sections, with airfoil type being NACA63-215.

The Wilson method is based on BEM theory. The method aims to maximize the utilization coefficient of wind energy. When designing the blade, it can take into account the influence of blade tip loss and airfoil lift-drag ratio on the optimum aerodynamic performance of the blade, so the designed blade has high wind energy conversion efficiency. Hence the incoming velocity is and the incoming angle can be expressed as follows:where and are the initial axial induction factor and angular induction factor, respectively. According to BEM theory, they can be deduced as follows:where is the blade solidity.

The acquisition of and can be iteratively repeated by the Wilson method as follows [11]:(a)Assume a reasonable set of initial values of and .(b)The incoming angle and attack angle α are calculated by (2a) and (2b).(c) and are obtained from the corresponding airfoil aerodynamic data library according to .(d)A new set of parameters of and is obtained by using (3a) and (3b).(e)The steps are repeated continuously until the two adjacent calculated values of and converge.

The final calculated values of and , as well as their curves fitted as expressions of general model Sin6 (Sum of six-order Sine) [15], are shown in Figure 2(a). Herein, stands for rotation radius of single cross-section.

According to the chord length formula [12], each chord length at distance , can be calculated aswhere is the number of blades. The curve can also be fitted as general model Sin6 expression (see Table 1).

At the same time, the linear density of the blade can be approximately estimated and fitted in Figure 2(b), according to the structure parameter and an expression of a density factor in terms of the air density in [15]. Also, the density factor has the corresponding fitting expression and fitting coefficients (see Table 1). In addition, considering the influence of the theories of blade cascades [16] on the attack angle, the optimal angle of attack can be obtained according to the Reynolds number. Then the pretwist angle is also obtained in Figure 2(b) according to (2a) and (2b).

The fitting curves of the 4 parameters in Figure 2 do not mean that they are more accurate than the optimized data itself. However, the equation of the fitting curve can be directly substituted into the aeroelastic system for integral operation. As mentioned above, the analysis and calculation process can be simplified as much as possible. Equation (5) is the expression of the fitted six-order Sin6 model of the related parameters, with the corresponding coefficients and the errors in the fitting process displayed in Table 1.where , , and are the corresponding structural items; for , , , , and ; for and .

After using graphical methods to evaluate the goodness of fitting, some items of the goodness-of-fit statistics for parametric models should be examined. The items include the sum of squares due to error (SSE), R-square, Adjusted R-square, and root mean squared error (RMSE), which can be found in Table 1.

It can be seen in Table 1 that the SSE values closer to zeros indicate that the models have smaller random error components and that the fitting will be more useful for prediction; the R-square values closer to 1 indicate that a greater proportion of variance is accounted for by the model, which means that the fitting explains greater proportion of the total variation in the data about the average. Also, the adjusted R statistics closer to 1 indicate a better fitting, with a RMSE value closer to 0 indicating a fitting that is more useful for prediction as an estimate of the standard deviation of the random component in the data.

2.3. The Fitted Aerodynamic Coefficients

Aerodynamics characteristics of airfoil are the basis of aerodynamic performance design of blades. Meanwhile the acquisition of airfoil aerodynamic data is often the most difficult part of blade design. The accuracy and completeness of the aerodynamic data obtained are of great significance to the analysis of the dynamic behavior of the blade. The blade of the wind turbine usually operates at the angle of attack range of −180°~180° with a wide range of Reynolds number. Hence the aerodynamic coefficients of airfoils in corresponding ranges are also required. Considering the actual situation of the present design, it is preliminarily estimated that the angle of attack will change between −90°~90°. The acquisition steps of aerodynamic coefficients are as follows [11]:(a)Firstly, the lift and drag coefficients within a range of small attack angles (−15°~15°) are obtained by using Xfoil software [11].(b)Using the worksheet “TableExtrap” in AirfoilPrep_v2p2 program that prepares the corresponding airfoil data for AeroDyn software [17], the range of small attack angles is extended to range of −90°~90°.(c)The validity of the fitting processes of structural parameters is demonstrated by the error data displayed in Table 1; the aerodynamic coefficients can also be fitted by the acquired data within the range of −90°~90° [15].

Figure 3 shows the curve of the aerodynamic lift coefficient (data ) versus the angle of attack (data ), as well as the aerodynamic drag coefficient (data ) against the angle of attack (data ), according to the acquired data from above steps (a~c). The structural items of coefficients ( and ) and the corresponding errors in the fitting process are also displayed in Table 1.

It should be emphasized that, in view of the actual situation of the design, only the aerodynamic coefficients were fitted to the angles of attack within the range of −90°~90°. If the large angles of attack and deep stall-induced situation are put into consideration, this design method can also be applied to the angles of attack within the range of −180°~180°.

2.4. Governing System and Solution

The 3D coupled blade body analysis method fully considers all the cross-sections along the spanwise length, and the cross-section displacements have complex coupling behaviors. Its solution needs strip theoretical decomposition, aerodynamic decomposition, and decoupling using Galerkin method [15]. The NCM method used here only extends the integral behavior of the 3D coupled blade body analysis method on the basis of the 2D typical section analysis method. Inserting (5) into (1a) and (1b) and considering the integral behavior along the blade spanwise corresponding to from 0 to R result in the equations governing the motions of the aeroelastic system as follows:where is the variable vector. Furthermore, in order to carry on subsequent analysis, (6) needs to be transformed into the state space expression as follows:where is the new variable vector.

Some cases intended to highlight the effects of the damping ratios and wind velocities on the vibration and stability of different sectional radius will be presented. Basic structure parameters are as follows:  rad/s,  rad/s, , . Figure 4 illustrates the real parts of both flap eigenvalues of vertical direction and lag eigenvalues of lateral direction at different positions (from 1/6R to R at intervals of 1/6R), under conditions of wind velocities from 5 m/s to 50 m/s at intervals of 5 m/s and damping ratios from 0 to 0.1 at intervals of 0.01. All the real parts of eigenvalues in both directions at every position are less than zeros, which means absolute stability. In fact, the eigenvalues here involve only homogeneous equation structure in the left-hand side terms in (7), which only concerns the mass matrix coefficients, damping matrix coefficients, and stiffness matrix coefficients in (1a) and (1b). If there is no aerodynamic force, most fundamentally, this implies constant stability of the structure. Therefore, in view of this characteristic, no matter how the change of wind speed and damping ratio, eigenvalues of homogeneous equation in (7) are also bound to be stable. Furthermore, in view of this characteristic, a novel pole assignment technology based on self-poles derived from these stable eigenvalues of homogeneous equation part in (7) is proposed, which is named as PPMO technology in Section 3.

As is known to all, the first step in the pole-placement design approach is to choose the locations of the desired closed-loop poles. The most frequently used approaches are to choose such poles based on experience and the quadratic optimal control approach. Otherwise, the system responses become very fast, with very large amplitudes; or the poles cannot balance the acceptable responses and the amount of control energy required. And the pole assignment technology based on self-poles derived from these stable eigenvalues can just overcome these shortcomings and save the trouble in finding the right poles.

3. Application of Pole Placement with Minimum-Order Observer (PPMO)

3.1. Pole Placement Based on Self-Poles

Different from specifying only dominant closed-loop poles in the conventional design approach [18], the present self-pole placement specifies all closed-loop poles including poles of two displacement variables and two velocity variables. The assignment of displacement poles can be efficiently performed owing to accurately measured displacements. However the assignment of velocity poles requires the inclusion of a state observer in the system. There is also a requirement for the closed-loop poles to be placed at arbitrarily chosen locations, which implies that the system should be completely state controllable. As depicted above, constant stability of the structure in (7) without action of aerodynamic forces decides the desired locations where the regulator poles can be placed, which is exactly what is said of the self-pole placement. The regulator poles are the eigenvalues of , where is determined by the control signal . The method of self-pole placement avoids the difficulty of selecting the desired poles and effectively improves the performance of the pole assignment, which shows the rapidity of the configuration process and the effectiveness of the configuration results. Present study shows that the selection of the desired closed-loop poles or the desired characteristic equation based on self-pole placement has not only the response speed superiority, but also the superiority in dealing with the sensitivity to disturbance and measurement noise.

According to Ackermann’s formula [19], the state-feedback gain matrix can be determined as follows:where satisfies the desired characteristic equation:

3.2. Minimum-Order Observer

During the pole-placement process, it is assumed that all state variables are available for feedback. In present study, velocity variables are unmeasurable for direct feedback and need to be observed. Notice that, in order to deal with velocity variables, one needs the derivative of displacement variables, which presents a difficulty due to the reason that differentiation amplifies noise. Hence it is important to estimate the velocity variables accurately. Because only two velocity variables need to be estimated, so the observer called minimum-order observer. Figure 5 shows the feedback system with a minimum-order observer.

The PPMO method exactly consists of the self-pole placement technology and the structure of minimum-order observer. Firstly, both the system and state variables of (7) are split into two parts, based on displacement variables and velocity variables , as follows:where are the state variables of the minimum-order observer system, then the feedback state variables of the transformation system in Figure 5 can be written as follows [19]:with the error () equation for the minimum-order observer described aswhere the state observer gain matrix is an 3 × 1 matrix and can also be deduced according to Ackermann’s formula as follows:where satisfies the characteristic equation and can be deduced as follows:where is the linear coefficient of corresponding characteristic equation.

The closed-loop poles of the PPMO system in Figure 5 comprise the closed-loop poles due to pole placement from the eigenvalues of matrix and the closed-loop poles due to the minimum-order observer from the eigenvalues of matrix (). Hence the pole-placement design and the design of the minimum-order observer are independent of each other. The characteristic equation of the PPMO system can be derived as . In other words, one step is to calculate the feedback gain matrix with a feedback law of , which has closed-loop poles at the specified values of all the self-poles; the other step is to compute a state-feedback matrix such that the eigenvalues of are those specified in the last half of the self-poles (corresponding to the two velocity state variables).

It should be emphasized that a necessary and sufficient condition for arbitrary pole placement is that the system be completely state controllable. In view of the constant stability of the homogeneous equation structure in (7), the system here is not only completely controllable which can be determined by the rank of the controllability matrix with , but also completely observable which can be determined by the rank of the observable matrix with .

4. Analysis and Discussions

4.1. Divergent Cases and PPMO Control

Figures 6(a)-6(b) show the uncontrolled displacements (a) and those controlled by PPMO (b) under the condition of wind speed of  m/s. The uncontrolled displacements exhibit unstable states, with vibration amplitudes more than the blade length of  m in a very small range of time, which is in fact an state of extremely divergent instability, from the point of view of numerical simulation. The controlled cases with PPMO control in Figure 6(b) present the states of convergent stability with very small vibration amplitudes, which demonstrates ideal control effects. Two characteristics of PPMO control performance need to be emphasized, one is small vibration amplitude as mentioned above, and the other is very small steady-state value which implies slight offset of vibration displacements of blade tip.

The linear quadratic regulator (LQR) is often used to analyze vibration control of rotating structure. In [20], the host structure with sensor-actuator pairs embedded along the spanwise direction was considered. The total output was fed to a controller and then two types of LQR methods were used to investigate rotating effects. Figure 6(c) illustrates comparisons of control effects between PPMO method and LQR algorithm. As can be seen, the LQR control results have relatively large offsets, which also means greater amount of control consumption which will be discussed later.

In order to verify that the PPMO control method can be commonly used at that given external wind velocities. The other cases of different wind speeds are investigated. Figure 7 shows comparisons of control effects between PPMO and LQR under conditions of wind speeds of  m/s (a) and  m/s (b), respectively. From the magnitude of the vibration, or from the steady-state value, almost the same conclusion can be obtained, which proves the robustness of PPMO control algorithm.

In addition, for the linear control strategy, the control cost (control consumption of energy) often affects the implementation of the control measures seriously. Therefore, the magnitude of the controller itself needs further exploration. Figure 8 shows amplitudes of PPMO controller and LQR controller versus different state variables, respectively. Herein, “ = 1, 2” refers to the first two displacement variables; “ = 3, 4” refers to the latter two velocity variables. It can be seen that the amplitude of PPMO controller is much smaller than the amplitude of LQR controller under various wind speeds.

From the point of view of the “control cost” of controller hardware mentioned above, another advantage of PPMO control over LQR control is that the implementation of LQR control requires full state-feedback, which will not only increase the number of the hardware sensors used for measurement, what is more, the measurement process itself of the state variables of velocities is a complex process.

In addition, the implementation of PPMO is an automatic process, including the configuration of self-poles which is also automatic. Performance of LQR control depends on the choice of weighting Q/R matrices which have no analytic solution. So the choice of Q/R matrices is artificial, which is another drawback of LQR control.

4.2. Feasibility Analysis of Hardware Implementation of Control Algorithm

In view of the large wind power system being mostly controlled by PLC, some intelligent control algorithms cannot be implemented in PLC hardware because of their complexity. At present, most of wind turbine control system algorithms employ fuzzy control, which can be implemented in PLC system. However, the fuzzy algorithm itself is tedious. For the other neural network algorithm, although the numerical simulation of control effect is more ideal, because of the complex integration problem, it is difficult to be implemented in the PLC hardware [21]. Hence the feasibility analysis of the hardware implementation of PPMO algorithm needs further discussion.

In present study, a hardware-in-the-loop simulation platform is built by OPC technology [21] between PLC controller system and MATLAB/SIMUILINK simulation environment to test the feasibility of hardware implementation. Figure 9 shows the experimental system, including the main interface of the human-machine interface (HMI), and the PPMO control interface, respectively. The experimental system consists of PLC controller system and MATLAB-OPC computer. PLC controller system runs the entire PPMO algorithm and sends the controller output to MATLAB-OPC computer to drive the aeroelastic system and accepts signal that is the output of aeroelastic system from MATLAB-OPC computer. The aeroelastic system model is run in MATLAB/SIMULINK environment. The signals of the control process can be displayed by HMI which is connected to PLC (CPU 224).

Figure 10 shows the touchscreen interface in HMI and displays the experimental signals of the vertical and lateral displacements without control (a) and those under PPMO control (b) based on the condition of wind speed of  m/s. Compared with the uncontrolled cases in Figure 6(a), Figure 10(a) shows almost the same change trends of vertical and lateral displacements without control. Not only can the trend of change be shown in HMI, but the specific values of change can also be shown, such as the lateral displacement in Figure 10(a). Also, Figure 10(b) demonstrates the controlled cases as depicted in Figure 6(b). Similar trends and magnitude values can be displayed, with similar conclusions obtained.

It should be noted that the data matching and parameter transmission planning of OPC technology are not discussed here, which is another specific and complex problem. What is more, the speed of PPMO control in experimental platform is faster, compared with the hardware implementation of conventional fuzzy control by OPC technology in [21], which might be due to the simplified matrix operations in PPMO algorithm.

4.3. Performance Improvement of Theoretical PPMO Control

Consider the controlled cases by PPMO in Figure 6(b), the vibration frequencies of both vertical and lateral displacements are relatively large. Especially for the vertical displacement in Figure 6(c), compared with LQR control, vertical vibration of PPMO exhibits greater frequency fluctuations, which is actually a real-time drawback. To improve the performance of PPMO control, an improved scheme to integrate PPMO method and LQR method together is highlighted with response effects demonstrated. The first step of the improved scheme is to determine the expected poles by the quadratic optimal control method [19], then to compare the expected poles with the self-poles of the system, and then to choose the best poles between them. The second step is to executive PPMO control programme on the basis of the best poles.

Still taking the cases of Figure 6 under wind speed of  m/s for example, Figure 11 demonstrates the vertical and lateral responses controlled by the improved scheme. The improved scheme, not only for vertical displacement but also for lateral displacement, not only relative to PPMO control, but also relative to LQR control, not only in the frequency fluctuations, but also in the steady-state values, even in the vibration amplitudes, has achieved greatly improved performance.

It should be stated that this performance improvement is only a theoretical improvement, and it only fully combines the advantages of PPMO control and LQR control and does not consider the shortcomings of LQR control. Therefore, this scheme has some unavoidable drawbacks in control cost, automation, and hardware implementation of control algorithms as mentioned above. Its implementation can only be expected in hardware performance improvement and control cost reduction in the future.

5. Conclusions

In this study, vibration control for vertical bending and lateral bending based on pole placement with minimum-order observer is investigated by numerical simulation and hardware-in-the-loop simulation platform. Some concluding remarks can be drawn from the results:(1)Structural modeling is investigated based on a NCM method between the 2D typical section analysis method and the 3D coupled blade body analysis method, which implies a compromise result that is bound to occur, but it can greatly reduce the difficulty of processing. Aeroelastic equation is realized by fitting of structural parameters and aerodynamic coefficients based on six-order Sin6 model, which can greatly simplify the integral behavior along the spanwise length and, at the same time, reduce the system variables (completely eliminate the aerodynamic variables).(2)PPMO method combines the advantages of pole placement and minimum-order observation, which is much better than LQR control either from vibration amplitude or from control cost. Pole assignment is based on self-poles, which is the most essential feature of PPMO control. Otherwise, the control performance of the proposed method is greatly reduced due to the difficulty of finding poles and the matching problem of configuration poles.(3)Feasibility analysis of hardware implementation of control algorithm is verified by hardware-in-the-loop simulation platform. The platform consists of PLC controller system connected with HMI, MATLAB simulation environment, and the OPC link between the two structures. In the experiment, the running speed of the platform and the accuracy of the displacement display further demonstrate the advantages of PPMO control.(4)In view of the fact that both simulation and experiment are the analysis and discussion of the control theory itself, some future work needs to be further envisaged. Practical engineering application of PPMO control can be realized via external pitch motion which can be fed back to the aeroelastic system. The pitch motion can be driven by a hydraulic system. Hydraulic systems have a rich variety of components and are easy to be controlled and should be easier to realize pole placement.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this article.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (Grant no. 51675315) and the Graduate Innovation Project of SDUST (Grant no. SDKDYC180333).