Abstract
Identification of a one-stage axial compressor system is addressed. In particular, we investigate the underlying dynamics of tip air injection and throttle activation to the overall compressor dynamics and the dynamics around the tip of the compressor blades. A proposed subspace system identification algorithm is used to extract three mathematical models: relating the tip air injection to the overall dynamics of the compressor and to the flow dynamics at the tip of the compressor blade and relating the movement of the throttle to the overall compressor dynamics. As the system identification relays on experimental data, concerns about the noise level and unmodeled system dynamics are addressed by experimenting with two model structures. The identification algorithm entails a heuristic optimization that allows for inspection of the results with respect to unmodeled system dynamics. The results of the proposed system identification algorithm show that the assumed model structure for the system identification algorithm takes on an important role in defining the coupling characteristics. A new measure for the flow state in the blade passage is proposed and used in characterizing the dynamics at the tip of the compressor blade, which allows for the inspection of the limits for the utilized actuation.
1. Introduction
Operation of axial compressor systems, as found in a number of aerospace applications, faces difficulties due to instabilities at the peak of its performance curves. To mitigate these instabilities, a number of different control mechanisms have been introduced. Among them are air injection into the tip of the blade passage [1–3] and the use of grooved casings [4–6]. Modeling of the overall compressor dynamics has been a challenging undertaking, resulting in the well-known Moore-Greitzer model. However, there does not exist a model relating the influence of air injection to the overall system dynamics. One of the aims of this work is to present a system identification based approach to characterize the influence of the air injection—commonly used to control the compressor—to the overall compressor dynamics.
As the system identification relays on experimental data, concerns about the noise level need to be addressed as well as unmodeled system dynamics. The unmodeled system dynamics can become an important issue when the experimental data do not include all system modes during the system identification experiment or the assumed model structure does not allow for accommodating one or more particular dynamic characteristics. We propose to utilize two different model structures to investigate the unmodeled system dynamics and the effects of the proposed optimization. These structures are the autoregressive moving average with exogenous input (ARMAX) model and the autoregressive with exogenous input (ARX) model. An innovation sequence is extracted by using an estimated parameter sequence. The ARMAX model and ARX model parameters are then utilized to construct a Hankel matrix, representing the sampled impulse response of the system relating tip air injection to pressure rise coefficient. A singular value decomposition of the Hankel matrix is performed for an eigensystem realization (ERA, [7, 8]) resulting in a state-space description of the coupling dynamics. However, prior to performing the state-space realization, a heuristic optimization algorithm is introduced to optimize the singular values with respect to the measured output data. This allows for minimization of the noise influence to the extracted state-space model as well as the effect the assumed model structure has on fitting to the recorded characteristics from the system identification experiment. The proposed system identification is used for data collected on an isolated rotor axial compressor system whose blade geometry allows for spike stall inception. The tip air injection is applied in a random fashion in order to allow for sufficient system mode excitation. The resulting pressure rise and associated flow coefficient are computed using a set of pressure sensors. The optimization of the Hankel matrix is accomplished by using a Tabu Search (TS) algorithm [9]. The TS searches the vicinity around each singular value defined by a percentage of the largest singular value. The eigensystem realization is done using a balanced realization, an input-normal form, and an output normal form realization in order to assess the impact of the optimized singular values. Another aim of this work is to characterize the dynamics within the blade passage at the tip of the rotor blade due to air injection at the leading edge of the rotor blade. Finally, the proposed hybrid system identification algorithm using TS is used to model the compressor dynamics excited by movements of the throttle and the resulting pressure rise changes.
In the following, the details of the system identification scheme as well as the proposed realization and optimization are introduced.
2. System Identification Algorithm
Consider a linear, time-invariant, discrete time system:where , , , and are the system matrices and , , are the output and input data sequences recorded during a system experiment with some sample frequency , resulting in discrete data points for each of the two sequences. is the number of outputs, and is the number of inputs; is the process noise and is the measurement noise, while is the state vector.
In order to estimate the state variable , a Kalman filter is used, provided the states are observable [10]. Hence, the system given by (1) and (2) can be reformulated. Defining the Kalman filter gain as where is the solution of the steady-state algebraic Riccati equation. Hence,whereThe system in (1) is in the process form, while the system in (4) is in the innovation form. Using (6) in (4) results and defining and , we arrive at the predictor form of the given system:The predictor form of (8) can be used to derive an AutoRegressive with eXogenous input (ARX) model: Note, is an asymptotically stable square matrix; hence if is sufficiently large, we haveEquation (9) becomesNote that are the system Markov parameters and are the Kalman filter Markov parameters. Also, the existence of is guaranteed if the system is detectable and is stabilizable [7].
Equation (11) is of the formwhere and . Assuming the system is observable and controllable, the following output vector can be created (for simplicity, we assume D = 0 or the first input = 0):Henceis the observability matrix. The system is observable if is of . To test for controllability, one can create a state vector:Hereis a square matrix. The discrete time system is controllable if and only if . Suppose the system of (4) and (5) is controllable and observable with rank n; then One notices that the Hankel matrix is composed of the system Markov parameters. Using the ARX model parameters , can be given asDefining the system Markov parameters asthen can be written asTo realize a state-space system of the form a minimum realization is utilized (see [7, 10]). The minimum realization allows for investigating the primary or dominant modes of the system to be identified.
Perform a SVD on , where and the singular values are ordered asTruncate and to have columns, which results into and . Hence the following equality exists:Hence, the Hankel singular values are the square roots of the eigenvalues of :Now there are different ways to interpret the realization. They are organized by what combination and association is being taken:(i) input-normal form:(ii) balanced form:(iii) output normal form: Choosing the balanced form, we haveHaving realized and , we can utilize their original definitionHence, is first columns of or of ; is first rows of or of , . To compute , one forms Expressing in terms of :Using SVD on the last equation yieldsFrom this, one can compute :
3. Estimation of ARMAX Model
Equation (12) is an ARX model. To incorporate a moving average term into (12), a third convolution term is added,where , , and are the model parameter matrices. The idea of using an ARMAX model is to accommodate the innovation and noise influence. The parameters of the ARX model can be estimated using standard least-squares techniques. For the ARMAX model, the following approach is adapted.
Definefor , where the subscript is used for denoting the time index. The innovation sequence can be estimated by using whereUsing and with , and definingwhere the subscript is used for the time index, the estimated ARMAX parameters are given byTo obtain a realization of the form of (39), the estimated model parameters of the ARX or ARMAX model are used to compute the system Markov parameters (see [11]) and construct a block Hankel matrix of the form given by (20) or (30). The realization then can be established by any of the three options given by (25), (26), or (27).
4. Tabu Search Optimization
The Hankel matrix in (20) is composed of the sampled impulse response of the system. In addition to the system’s response, there is noise, unmodeled system dynamics, and responses from auxiliary systems that are coupled with the system to be identified. As there is no method to separate the given impulse response from coupled noise, unmodeled dynamics, and auxiliary system responses, we propose to imbed an optimization into the eigensystem realization (ERA) algorithm given by [8]. The optimization is done using an Enhanced Tabu Search (ETS) algorithm. The ETS employs a two-stage process, where during the first stage a list of promising areas in the search field is established. During the second stage, the promising areas are further searched by a regular TS algorithm. Details on the employed ETS and TS are found in [9]. The ERA is updated such that after selecting the number of states of the system—based on the magnitude of the singular values in —the entries of the selected singular values are searched within a close vicinity of their nominal value established by the SVD. The search areas are defined as a percentage of the largest singular value. The cost function for the ETS and TS is given bywhere is the estimated output from the resulting realized system. An easy extension to the given cost function is to incorporate a priori information of the system. For example, if the fundamental natural frequency of the system to be identified is known, (40) can be extended to where are weighting factors and is the resulting estimated natural frequency.
5. Experimental Setup
In this work, we utilize a one-stage compressor system with a blade geometry that allows for spike inception. The proposed identification scheme is used to identify the dynamics at different flow conditions, including near stall conditions. The dynamics of the compressor is given by the changes within the blade passage area as well as the changes in the pressure rise coefficient computed by the input and output pressures. Hence, we extract three separate models, one for the dynamics within the blade passage due to air injection, one for the input/output relationship of the compressor due to throttle movements, and one that characterizes the dynamics due to air injection pressure as the input of the system and the pressure rise coefficient as the output of the system. The compressor employed for this research is a low-speed rotor with characteristic parameters as given in Table 1. The compressor exhibits stall at 40–50% of its rotating frequency. The operating range of this compressor is given by a flow coefficient between 0.58 and 0.49, where the latter bound represents the vicinity of the stall point. The experiments were carried out using a smooth casing; that is, there are no grooves within the tip clearance casings. The average tip clearance for the setup is 1 mm. A number of pressure sensors are used, as shown in Figure 1.
Sensors 1 and 8 are utilized to compute the flow coefficient and pressure rise. A separate sensor is mounted at the injection port to measure the injection pressure. The pressure sensors measure at a sampling rate of 20,000 Hz. There are eight air injectors, equally distributed over the circumference, injecting air at an angle of 15 degrees with a steady pressure. Sensors 2 to 7 are dedicated to the dynamics within the blade passage.
For the actuation, there are two inputs. One input is given by the movement of the throttle valve, which is used to set the operating point. The other input is accomplished via high-pressure air injection using the eight injectors equally distributed around the circumference of the compressor annulus. The setup of the measurement and injectors is shown in Figure 2(a) while the throttle valve cone movement setup is shown in Figure 2(b).
(a)
(b)
6. Results and Discussion
6.1. Identification Results for Overall Dynamics due to Injection
We first present the results for inferring the compressor dynamics modeling the dynamics related to the air injection and overall output of the compressor as measured by flow coefficient and pressure rise coefficient. Considering Figure 3, the system identification results for a flow coefficient of 0.55 are shown by the extracted Bode plots of the various realized systems. The original ARX and ARMAX model, using a balanced realization and no optimization of their singular values, indicates that both have captured very similar characteristics.
The first fundamental frequency is at around 557 Hz. Once optimization is used as well as different realization forms (input (IN), output (ON), or balanced (BN) normal), the Bode plots of the two models (ARX and ARMAX) start to diverge from each other. For the ARX based model set, the frequency responses converge to a very defined Bode plot, with a fundamental frequency of 477 Hz. For the ARMAX based model set, the resulting frequency responses are reduced in their emphasis to the fundamental frequency. The ARX based model set also shows a much improved magnitude for the first fundamental frequency. Regardless of the realization (balanced, input, or output normal), all of the optimized ARX based models seem to agree with the first fundamental frequency and its amplitude. The phase plot does not show much change from unoptimized to optimized system identification results, for both model structures and all realizations.
At a flow coefficient of 0.55, the compressor operates sufficiently far away from the stall inception. The optimization of the singular values for the system identification algorithm is based on minimizing the error between the simulated model output and the measured overall system output, as given by (40). To gain some understanding of what these system identification results mean, we shall consider the difference of the two model structures, that is, ARX and ARMAX. The ARX model can be given in system block formulation as shown in Figure 4.
The corresponding system equation is given bywhere is the backshift operator and and are polynomials of the numerator and denominator. The ARMAX model is given in Figure 5.
The corresponding system equation is given bywhere is the corresponding polynomial responsible for modulating the noise sequence . Comparing (42) with (43), we notice that the ARMAX model has a built-in capacity to model the unobservable noise term . However, in this case, we not only have an unobservable noise term, but also the overall system dynamics, apart from the injection induced dynamics. When using the TS optimization, the incorporated cost function entails the overall output, that is, the system dynamics output, unobservable dynamics, and the dynamics resulting from the noise sequence. As the system identification algorithm has no means to separate the different dynamics and the noise influence, it models all of it with the help of the Moving Average (MA) portion of the model structure. This is corresponding to the polynomial in (43). As the input/output energy of the collected data is the same for all experiments, the modeling of the MA portion drains some of this energy away from the coupling dynamics and results into lower magnitude plots in the Bode diagram for the ARMAX models. The measured data is filtered prior to the application of the proposed optimization, using a notch filter to reduce the influence of the blade passing frequency (BPF). The sampling time of the data collection is set to 20 kHz, making the resulting Bode plot span over a large bandwidth. However, the overall dynamics of the compressor is found to be at much lower frequencies. This is given by the stall frequency of this compressor to be 17 Hz. Hence, the discussed frequencies in the Bode plot likely do not represent the coupling dynamics. From measurements at various locations through the compressor ducting, it was noticed that the 477 Hz frequency exists everywhere (sensors 1 and 8 in Figure 1). The true meaning of this frequency is not understood at this time. However, observations from experiments indicate that the amplitude of this frequency decreases as the throttling increases. This observation fits the characteristics of the two Bode plots given by Figures 3 and 6.
Without any current theory explaining the observation of the mode at 477 Hz, our focus will be directed towards the low frequency behavior of the system. Doing so, we assume that the coupling dynamics are found closer to the slow dynamics frequency region. Utilizing air injection as a control input breaks up the flow structure at the blade level. The overall dynamics is given by the flow coefficient and pressure rise coefficient. These two measures seem to reside in the low frequency area of the data. Hence, a second filter is proposed to be used to prepare the collected data. In particular, a low pass filter with a cutoff frequency of 150 Hz is utilized. The filter employs 20 filter weights and a Kaiser window. In addition, the data is downsampled from 20 kHz to 400 Hz. The filtering and downsampling are used to ensure that the system identification and its optimization focus on frequencies below 150 Hz. The proposed optimal identification algorithm is employed to extract dynamic relationships between the input (injection) and the output (flow coefficient). For system identification, the data is also prepared by subtracting the dc offset; that is, the data has a zero mean. For the results, this implies that the output of the models expresses the change in flow coefficient, rather than the flow coefficient itself.
In the following, Bode plots of the various identified models are given, using the proposed filtering and resampling. Figure 7 depicts the models for a flow coefficient of 0.58 (far away from stall). From Figure 7, one can recognize two peak frequencies for the ARMAX model based systems. The first fundamental frequency is at 79.6 Hz; the second is at 114.6 Hz. For the ARX model based systems, the first fundamental frequency is at 81 Hz; the second is at 125.7 Hz. Figure 8 depicts the results for a flow coefficient of 0.55. The identification resulted in some unstable models for the ARMAX based systems; hence the results for these (blue lines) are disregarded for this flow coefficient. For the ARX based systems, the first fundamental frequency is found at 77.9 Hz, and the second at 122.5 Hz. Finally, Figure 9 depicts the outcome of the identification experiments for a flow coefficient of 0.51 (close to stall). For the ARMAX based models, the first fundamental frequency is found at 79.6 Hz and the second at 113 Hz. It is interesting to note that, for the ARX based models, three fundamental frequencies are found. In particular, for the ARX based models, the first fundamental frequency is found at 44.6 Hz, the second frequency at 81 Hz, and the third at 114 Hz. These values seem to be close to the rotor frequency and its harmonics. However, the rotor frequency and BPF are filtered out using notch filters. At a flow coefficient of 0.51, the simulated output from the identified model and the actual output (change of flow coefficient) are plotted in Figure 10.
A test to see if the extracted models can predict the output should provide indication if these frequencies are based on residuals of the rotor frequency or if they actually model the coupling dynamics between the injection and the pressure rise coefficient. Comparing Figures 10 and 11, the effect of the optimization is shown by the smoother output of the simulated model due to optimizing the singular values, and hence reducing the noise influence.
A similar observation can be made for the ARMAX models shown in Figures 12 and 13. Here, the ARMAX model using no optimization for the realization shows marginal stability properties (Figure 13). The performance of the ARMAX model is much improved and yields a stable system, shown in Figure 12. When overlaying the injection pressure measurements, as shown in Figure 14, the correlation between the injection pressure and the ARX model outputs is evident. The ARMAX model does capture this behavior too; however the output is inverted. Comparing the simulated output (Figure 10) with the injection pressure (Figure 14) one could conclude that the model picks up the input as a direct transmission term, representing the coupling in a static fashion. Noticing the scale of the simulated and measured response, the extracted dynamics is about half of the measured output in terms of magnitude. The optimization of the singular values seems to improve the performance of the extracted models (i.e., more smooth and clean signals):
The system given by the optimized ARX model (balanced realization) can be given in state-space form (see (1) and (2)) as shown above. Looking at the term, the direct transmission term, the magnitude is comparable to the elements in the other state-space matrices. Hence the system is not static; that is, the output is a combination of the dynamic part (the , , and matrices) and the direct transmission matrix .
The ARX/ARMAX model order is set to 12. As these singular values are optimized using a TS algorithm, the magnitude of the cost function—as defined by (40)—is plotted versus the number of iterations. This is shown in Figure 15.
For the simulations, a total of 45 Markov parameters were used to construct the block Hankel matrix given in (20) and (30). The computation of the ARMAX model used , and the TS algorithm is run for 200 iterations. A Tabu list and promising list of length 10 is used, that is, allowing 10 regions to be included in the intensification portion of the ETS algorithm, while keeping the last 10 steps of the search as forbidden moves.
6.2. Identification Results for Throttle Movement and Overall Dynamics
The proposed optimization embedded in the identification algorithm and the different cost functions are first tested using simulations. A simple second-order system with a natural frequency of 17 Hz is utilized. The simulations are carried out using different process and measurement noise levels. For each case, 20 simulations are used to statistically characterize and compare the performance of the proposed identification routines. The ETS algorithm uses 200 iterations for each simulation; the model order of the ARX and ARMAX models is . For the construction of the Hankel matrix (30), 75 Markov parameters are used. The search area around the singular values is set to be 10% of the nominal value of the largest singular value. The computation of the residual sequence for the ARMAX model uses two steps; that is, . Table 2 lists the results of the simulation using noise values of 1% and 5% standard deviation. For each approach, the nominal value (nom.), the optimized value using (40) (opt.), and (41) (opt. ) are provided as the mean value and the standard deviation. For both cost functions, the weighting coefficient is set as .
The error squared and correlation coefficients are computed based on 1000 validation data points. In Table 2, the two best performing algorithms’ results are highlighted. It is evident from Table 2 that the error of the ARMAX model is greatly reduced using the proposed optimization. The results of the ARMAX model and ARX using prior knowledge of the natural frequency indicate having more consistent results based on their lower standard deviation values compared to the standard ARX model. Note also that the error and correlation incorporate the noise as part of the signal. Hence, using the methods which do not incorporate the knowledge of the natural frequency may be prone to model the noise by overfitting. Using a higher noise level, Table 3 presents the simulation results with 10% and 15% noise level.
From the inspection of these results, we conclude that the ARX and ARMAX based model identification with optimization either with or without prior knowledge of the natural frequency works sufficiently well for use of identifying the overall system dynamics of the compressor system. Considering the overall dynamics of the compressor system between the input (throttle movements) and the corresponding pressure rise coefficient, the proposed system identification algorithm with prior information of the natural frequency is used. The identified model fits well for the three different flow coefficients of 0.51, 0.55, and 0.58, indicating that the overall dynamics does not change noticeably during change of operating point. It is expected that this may change, as the operating point gets closer to stall. For the identification, the same parameters are used as given in the simulation section. The resulting natural frequency of the extracted model is 17.1 Hz using the ARMAX base model with optimum and 16.8 Hz for the ARX base model with optimum , that is, using (41). The identified model in discrete time is given in the Appendix. Comparing the results of the nonoptimum methods with the resulting identified model, the natural frequencies are well off the known frequency of 17 Hz. For example, the ARX based system identification yields a natural frequency of 2,191 Hz. The incorporation of the a priori information also helped the correlation of the estimated output to improve considerably compared to the nonoptimum identification method treated in this work.
6.3. Identification Results for Injection and Blade Tip Flow Dynamics
For the identification of the dynamics within the blade passage, sensors 2 to 7 are used to capture the resulting dynamics at the blade tip area. As the sensors measure the pressures at all times, the data does not correspond to a single blade passage (the compressor runs at 2400 rpm). A hall effect sensor is used to phase lock the data and computes the pressure distribution within the blade passage. A sample pressure distribution is given in Figure 16 for a flow coefficient of 0.58. The computation of the pressure distribution as shown in Figure 16 requires averaging signals over hundreds of rotations.
For the system identification of the dynamics within the blade passage, the phase locking approach will not provide sufficient information, as only approximately 2% of the data can be associated with one rotation of the blade passage. We assume that the measured data of one blade passage corresponds to the state of a set of blade passages at that time. This is not to say that all blade passages have the same flow distribution at a given time. We rather state that the condition of the flow is similar. By using this assumption, we propose to use an information based entropy measure for the state of the dynamics within the blade passage. Recent work utilizes the correlation coefficient between the pressure data of the blade passage and the pressure data corresponding to the same blade passage one revolution prior. The correlation coefficient is used to characterize the state of flow in a blade passage. The lower the correlation coefficient, the more unsteady the flow. By using an entropy measure, there is no need to utilize data from over one rotation of the compressor. This entropy can be computed directly with the data present and hence has the same sampling frequency as the input. Therefore, we can associate the state of fluid flow with one characteristic number. The Shannon entropy is commonly used in information theory and measures the amount of information contained in a message. Here, one uses the entropy of a signal and classifies its predictability by a low entropy value and its randomness or amount of disorder by a large entropy value. The entropy is computed aswhere is the number of bins for the construction of the histogram, is the set of data to be used for the computation, and is the bin probability. An extension to the Shannon entropy is the spectral entropy, which can be used to characterize the distribution of energy in a signal. The spectral entropy can be computed aswhere is the Power Spectral Density, is the frequency, and is the normalized PSD.
By assuming equivalent flow condition at a given time for a set of blade passages, we artificially increase the spectral entropy value, due to the variation of flow energy among different blade passages at a given time. The entropy can be assessed at each sensor location and hence provides a snapshot of the flow characteristics. This is shown schematically in Figures 17 and 18 which depict the computed spectral entropy of the leading sensor (Sensor 2) using (46). For the first part of the time history, no injection is used, while, for the second part, injection is used in some random fashion (on/off). In information theory, often the zero frequency component of the spectral entropy is disregarded in order to measure a signal’s true disorder. Figure 18 depicts the spectral entropy with (red) and without (blue) the zero frequency component. For the purpose of system identification, from this figure, we can safely assume that the inclusion of the zero frequency component has no effect on the inferred system dynamics. It should be noted that the disadvantage of using an entropy measure as the measured quantity—for characterizing the flow—is that it increases the computational cost compared to the correlation coefficient approach and hence is only useful for offline work.
Figure 19 shows the spatial distribution of the spectral entropy within a blade passage using the seven sensors as given in Figure 2. The entropy values are given for two cases, without air injection and with air injection, for three different flow coefficients. From Figure 19, it appears that the injection affects primarily the entry of the blade passage where the spectral entropy level is raised, primarily for a flow coefficient close to stall (flow coefficient = 0.51). There is little effect due to injection into the end of the blade passage for all other flow coefficients.
Utilizing the optimized identification approaches for the ERA with ARX and ARMAX base models, the identification results for the dynamics of the flow at the entry of the blade passage are reduced to the fundamental natural frequency and given in Table 4.
The identification utilized 45 Markov parameters to construct the Hankel matrix given in (30) resulting in state-space systems with three states, for each flow coefficient. The optimization parameters are the same as those given for the simulation work stated earlier. Considering the resulting natural frequencies of the extracted state-space models, it is evident that these frequencies greatly reduce in magnitude as they approach the stall limits. The stall frequency is between 40% and 50% of the rotating frequency (~40 Hz). As future work, the convergence to this stall frequency will be investigated.
7. Conclusions
In this work, we propose a system identification algorithm utilizing a TS based realization. The TS optimization along with the system identification algorithm is used to investigate two different model structures and their characteristics applied to an axial compressor system. In particular, three dynamical models are extracted which characterize input-output behavior at different levels of the compressor system. One of these models describes the relationship between air injection at the blade tip and the resulting overall compressor dynamics represented by the computed change in flow coefficient. Considering the presented results: injection at the tip of the leading edge of the compressor blade has an influence on the overall dynamics of the compressor. These types of coupling dynamics may be of interest in developing more efficient control schemes for compressor control. Such controllers can contribute to extending the stall margin improvement (SMI) and run axial compressors at higher efficiencies. Another model extracted using the proposed identification algorithm captures the dynamics between the throttle movement and the overall dynamics as measured by the change in the computed flow coefficient—using measured data points. The proposed identification scheme embeds prior knowledge of the compressor and therefore guides the identification to more accurate results. Finally, the proposed identification algorithm is used to infer a model relating the air injection—treated as an input—at the blade time to the resulting flow dynamics near the tip of the compressor blade. At this stage, such models relating the air injection to flow characterizations near the tip may help in gaining a better understanding of the effect from the air injection, its reach into the blade passage, and its potential stall dynamics at this location. The proposed identification scheme with an embedded TS optimization algorithm as presented in this work seems to be able to aid in developing any of the three models stated above.
Appendix
Discrete time state-space model characterizing the dynamics between throttle movements and pressure rise coefficient changes is
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was supported by a generous grant from the National Science Foundation of China on Project no. 51306178.