#### Abstract

With depletion of traditional energy and increasing environmental problems, wind energy, as an alternative renewable energy, has drawn more and more attention internationally. Meanwhile, wind is plentiful, clean, and environmentally friendly; moreover, its speed is a very important piece of information needed in the operations and planning of the wind power system. Therefore, choosing an effective forecasting model with good performance plays a quite significant role in wind power system. A hybrid CS-EEMD-FNN model is firstly proposed in this paper for multistep ahead prediction of wind speed, in which EEMD is employed as a data-cleaning method that aims to remove the high frequency noise embedded in the wind speed series. CS optimization algorithm is used to select the best parameters in the FNN model. In order to evaluate the effectiveness and performance of the proposed hybrid model, three other short-term wind speed forecasting models, namely, FNN model, EEMD-FNN model, and CS-FNN model, are carried out to forecast wind speed using data measured at a typical site in Shandong wind farm, China, over three seasons in 2011. Experimental results demonstrate that the developed hybrid CS-EEMD-FNN model outperforms other models with more accuracy, which is suitable to wind speed forecasting in this area.

#### 1. Introduction

With the rapid development of China’s economy over the last few decades, China is becoming a big power in energy consumption. The consumption of coal, crude oil, and natural gas constituted 68%, 19%, and 4.4% of the total annual primary energy consumption in 2010, respectively, and the corresponding proportions have changed to 66.6%, 18.8%, and 5.2% in 2012 [1]. The total annual primary energy consumption has increased from 3249.39 million tons of standard coal as calorific value calculation in 2010 to 3617.32 million tons in 2012 [2]. Fossil fuels are nonrenewable resources; moreover, a tremendous amount of greenhouse gases is exhausted, which inevitably increases the burden of global climate warming [3]. Burning fossil fuels is responsible for environmental issues that are high on the political agenda these days. Thus, in order to deal with the energy crisis and global climate warming, clean and renewable energy is urgently needed.

With depletion of traditional energy and increasing economic development, the ecological environmental problems have become increasingly acute in recent years. And many countries in the world are struggling to look for new alternative energy sources. The potential role of renewable energy sources, such as wind, solar, geothermal, tidal, biomass, and hydropower, is becoming increasingly significant and has received increasing attention as they offer various advantages over nonrenewable, conventional energy sources in terms of environmental health and safety [4, 5]. Wind energy, as an alternative renewable energy to fossil fuels, has drawn more and more attention, which is plentiful, widely distributed, clean, and environmentally friendly because of producing no greenhouse gas emissions during operation and exploitation. With large scale development and utilization, wind energy has become one of the most rapidly growing renewable energy sources and has been regarded as an appealing alternative to conventional energy generating from fossil fuels [6]. The effects on the environment of wind energy are generally less problematic than those from other power sources; thus wind energy is very important among the low-carbon energy strategies, which has the potential to achieve sustainable energy supply.

Over the last decade, wind energy has developed very fast worldwide. At the end of 2009, the global wind energy’s installed capacity reached 158 GW (gigawatts). In China, as a big energy consumer, the wind power generation capacity has been growing rapidly. And in 2010, China replaced the United States as the biggest wind market. The current total capacity of wind farms in China is 91.413 GW in 2013 [7]. The integration of wind energy provides a number of advantages, but the large scale penetration and high-efficiency utilization of wind power pose many challenges as well in power system operations and planning, primarily because of the uncertainty and intermittency of wind speed [6]. The wind speed is random and intermittent, resulting in that wind power has a strong volatility and uncontrollability. In the electricity system the power supply must be in correspondence with the power demand all the time. However, the variation of wind power output makes it difficult to maintain this balance. To use the increasing wind energy effectively and safely, one of the possible solutions to the balance challenge is to improve the wind speed forecasting precision [8, 9].

In the last several decades, different wind speed predicted models have been developed to realize wind speed forecasting including linear time series analysis approaches ARIMA models [10, 11], Kalman filters [12, 13], and machine learning techniques, involving support vector machines [14, 15] and a variety of artificial neural networks such as multilayer feedforward neural networks [16–18] and recurrent neural networks [19, 20]. Despite signs of some improvements in these methods, however, there is always some limitation that each type of model directly tries to capture only linear patterns or nonlinear patterns. In order to get the most of each method, some hybrid models emerge as the times require which mostly combines two or more of the above models [21–23].

Zhang et al. [24] employed four improved adaptive coefficient approaches optimized by particle swarm optimization (PSO) to forecast the daily mean wind speed, of which the simulated results showed the PSO obtained an observable improvement in forecasting performance. The fractional-ARIMA models [11] were applied to forecast hourly average wind speed on the day-ahead (24 h) and two-day-ahead (48 h) horizons, which parsimoniously capture the presence of correlations of the wind speed. Then, the forecasted wind speeds were used in conjunction with the power curve of an operating wind turbine generator to achieve corresponding forecasting output powers. A hybrid ARMA-GARCH method [25] was presented to forecast hourly wind speed based on a series of 7-year data in Colorado, USA, the results of which showed the performance of the ARMA-GARCH was satisfactory. Guo et al. [26] proposed a modified EMD-FNN model combining empirical mode decomposition (EMD) with feedforward neural network (FNN) ensemble learning paradigm which performed better accuracy than that of the basic FNN and unmodified EMD-FNN. Cadenas and Rivera [27] developed a hybrid ARIMA-ANN model, in which the ARIMA was used to forecast the wind speed, and then in order to reduce the final errors, ANN was built to capture the nonlinear tendencies of the obtained errors of ARIMA model. An ARIMA-ANN model and an ARIMA-Kalman model were compared by Liu et al. [28] for multistep ahead prediction of wind speed, in which the ARIMA was not employed to predict the wind speed directly but was adopted to choose the best parameters for the ANN and Kalman components. A hybrid approach [29] based on the ensemble empirical mode decomposition (EEMD) and the support vector machine (SVM) was proposed to forecast the mean hourly wind speed, which proved to be effective in improving the prediction accuracy. Guo et al. [30] proposed a hybrid wind speed forecasting method based on BP neural network and eliminating seasonal effects using seasonal exponential adjustment, which forecasted the daily average wind speed more accurately than the BP model without adjustment. Another novel hybrid approach [31] combining seasonal adjustment method (SAM), exponential smoothing method (ESM), and radial basis function neural network (RBFN) was employed to improve the quality of wind speed forecasting.

The rest of this paper is organized as follows. In Section 2, the development of wind energy in China is described as well as a brief analysis to the available wind speed data in Shandong wind farm, China. The novel hybrid CS-EEMD-FNN method is presented in Section 3, which combines fuzzy neural network (FNN) with ensemble empirical mode decomposition (EEMD) and cuckoo search (CS). Statistics measures to determine the accuracy of the forecast are presented in Section 4. Section 5 gives the final forecasting results in the studied area. The performance of four methods is compared with the under four error evaluation criteria in Section 4, and the best model selected in Section 5 is applied to the wind speed prediction in the practical application of wind farms. Finally, Section 6 concludes the findings of this paper.

#### 2. Wind Energy in China

Wind energy has developed very fast worldwide over the last decades, especially in China. The new global total installed capacity of wind power at the end of 2013 was 318.105 GW, representing cumulative market growth of more than 12.5%, an excellent industry growth rate given the economic climate, even though it is lower than the annual average growth rate over the last 10 years of about 21%. China, the largest overall market for wind since 2009, had a good year and once again gained the top spot in 2013 with an installed capacity of 91,412 MW, sharing 28.7% of the global total.

##### 2.1. Development of Wind Energy in China

With large land mass and long coastline, China has exceptional wind resources. In 2007, the total installed wind power capacity in China (excluding Hong Kong, Macao, and Taiwan) was only 5.85 GW, and since then, the scale of China’s wind power industry has continued explosive growth. In the year of 2008, China’s new wind power installed capacity reached 6.15 GW, achieving a new capacity growth rate of 108% and the cumulative installed capacity jumped over 12 GW [32]. In 2011, the new annual installed wind power capacity became 17.63 GW, and by the end of 2011 the cumulative installed capacity was 62.36 GW, almost 11 times that in 2007. In 2011, China was the world’s second-largest wind producer, generating 73.2 billion kWh (kilowatt-hours), a level about 64% higher than in 2010 when the wind-generated electricity was 44.6 kWh [33]. In 2012, wind-generated electricity in China amounted to 100.4 billion kWh, accounting for 2% of the country’s total electricity output, up from 1.5% in 2011 [34]. Wind power generated 134.9 billion kWh of electricity in 2013, contributing 2.6 percent of the country’s total electricity generation. The Chinese wind capacity rose from 75.3 GW in 2012 to reach 91.4 GW by the end of 2013 with an increase of 21.4%, cementing the global lead in terms of cumulative installed wind power capacity. All the statistic data, including “Cumulative capacity” and “Annual installed capacity” and the corresponding growth rates, are shown clearly in Figure 1. In China, Inner Mongolia, Xinjiang, Liaoning, Shandong, and Guangdong are rich in wind resources and the wind power industry has developed very rapidly.

##### 2.2. The Available Wind Speed Data

As located in demand centers of China, Shandong wind farms have got good opportunities for harnessing the power of wind; meanwhile, the significant wind power potentials of the wind farms are around the coastal areas and islands. The total installed capacity of wind power in Shandong Province is 4562.3 MW in 2011, which becomes 6980.5 MW in 2013 ranking the fourth in China. Moreover, Chinese government has formulated a series of preferential policies for the development of wind electricity in this area; thus, accurate wind speed forecasting is of crucial importance for planning and operating a power system of the city’s sustainable development. So case studies are carried out using wind speed data from Shandong wind farms of China.

The wind speed data used in this paper are collected from Wind Turbine Generator of number 25 to number 28 (WTG25–WTG28) in Shandong wind farm of China. The studied time range covers January 1, 2011, to July 31, 2011. The ambient wind speed data are sampled with a certain time interval of 10 minutes, so there are 144 data records in a day. Figure 2 shows the variation trend of wind speed throughout the studied time range. Then we implement multistep forecasting for the 10-minute wind speed using the above-mentioned four models (FNN, EEMD-FNN, CS-FNN, and CS-EEMD-FNN), respectively.

#### 3. Methodology

##### 3.1. T-S Fuzzy Neural Network (FNN)

T-S fuzzy system [35] was initially proposed by Takagi and Sugeno, which could be presented as the linear combination of input variables. T-S fuzzy model expresses the knowledge with fuzzy and qualitative features and thus can deal with the uncertainty, nonlinear, and ill-posed problems. When T-S fuzzy system is applied into general fuzzy neural network, it brings T-S fuzzy neural network. Recent results [36, 37] show that T-S fuzzy neural network is a promising approach to capture the advantages of both fuzzy systems and neural networks and discard their weakness. Moreover, a fuzzy neural network can self-adjust the parameters of fuzzy rules using neural network based learning algorithms. Now it has been widely and successfully used in a variety of applications in many fields [38–40].

###### 3.1.1. The Common Structure of T-S Fuzzy Neural Network

T-S fuzzy system consists of a set of fuzzy IF-THEN rules. The rules of T-S fuzzy system can be represented as follows.

Rule : If is , is , and is ,Here () is the th fuzzy set of input in the fuzzy system, () is parameter of the fuzzy system, and is the output through the fuzzy rules, which is expressed as the linear combination of input variables.

The architecture of the T-S fuzzy neural network (FNN) is a four-layer fuzzy neural network shown in Figure 3. The four layers are the input layer, the fuzzification layer, the fuzzy rules calculation layer, and the output layer. Function of every layer is described as follows.

*Layer 1: The Input Layer.* The input layer is connected with the input variables. Let denote the input variables. The node number of this layer is . Each neuron in this layer represents an input variable , .

*Layer 2: The Fuzzification Layer*. Each neuron in this layer represents an If-part (or premise part) of a fuzzy rule. The function of the second layer is to calculate the fuzzy membership function of each input variable. The corresponding membership function of the input subvariable is calculated aswhere is the th membership function in the th neuron; is the center of the th membership function in the th neuron; is the width of the th membership function in the th neuron; is the number of input variables; is the number of the fuzzy subsets (the fuzzy division number of ).

*Layer 3: The Fuzzy Rules Calculation Layer*. The third layer calculates the applicability of each rule, which computes membership grade according to each fuzzy membership function. The degree that the input variable matches rule is computed by products of the grades of membership functions:where is also called the firing strength of rule .

*Layer 4: The Output Layer*. The fourth layer takes on the calculation of the output of a neuron converting the results into unitary form. Each neuron in this layer represents an output variable as the weighted summation of incoming signals. The output of a neuron in the layer 4 is

###### 3.1.2. The Learning Method of T-S Fuzzy Neural Network

(a) The error cost function can be defined asHere and are the expected output and real output, respectively. In FNN the gradient method is used for optimization and the Gauss function is chosen as the fuzzy membership function; then the learning parameters can be got by partial derivative method. These parameters are the weights of fuzzy system , the center values , and the width values of the membership function of the nodes in the second layer of FNN.

(b) The learning recursive formula of can be represented asHere , , is the rate of learning. is the number of sampling:When fixing the parameter , the whole structure is similar to the network of the typical fuzzy BP network essentially.

(c) The learning formula of and can be described as

##### 3.2. Ensemble Empirical Mode Decomposition Method (EEMD)

The ensemble empirical mode decomposition (EEMD) method is improved in order to overcome the shortcoming of the mode mixing problem in the empirical mode decomposition (EMD) method, which was accepted as a method of data preprocessing and was initially developed in 1998 [41]. Nowadays, the EEMD method has been successfully employed to remove the noise embedded in the time series in nature science and social science researches [42–44].

The EEMD has been widely used to decompose nonlinear and nonstationary time series data into a finite and small number of intrinsic mode functions (IMFs) through the sifting process using the Hilbert-Huang transform (HHT) until the final data series are stationary. An IMF can be best defined as a hidden oscillation mode that is embedded in the data series, since it is allowed to be nonstationary and be amplitude or frequency modulation [45]. Thus, the EEMD approach is a significant and meaningful part of data cleaning to the operation of forecasting wind speed, in which the highest frequency oscillation mode is identified as the first IMF, and the frequency of the following IMFs decreases on each step. All the IMFs should satisfy the following two conditions: (a) the number of extremes (sum of maxima and minima) and the number of zero-crossings must be equal or differ at most by one in the whole data set; (b) at any point, the mean value of both the local maxima envelope and local minima envelope is zero. With the above definition of IMFs in the EMD method, any original signal can be decomposed as follows.

Identify all the local extremes (maxima and minima) of , and then connect all the local maxima by a cubic spline line to form the upper envelope of , and similarly, the lower envelope is produced with all the local minima by repeating the procedure. The upper and lower envelopes should cover all the data of . Their mean is designated as where , and the difference between the signal and is the first component ; that is,

Theoretically, should be an IMF satisfying the above two conditions of an IMF. However, that is usually not the truth because changing a local zero from a rectangular to a curvilinear coordinate system may bring about a new extrema; thus further adjustments are needed [26]. Therefore, the above sifting procedure must be repeated again.

In the subsequent sifting processes, can be treated as the data needed to be decomposed in the next iteration:

The sifting process will be repeated times, until becomes a true IMF satisfying the above two conditions; that is,where is the mean of the upper envelope and the lower envelope of and . Then, it is denoted by

Moreover, a standard deviation between the two successive sifting operations is defined to determine whether the sifting process should be stopped, which is expressed as

Particularly, can be accepted as the agreed termination condition, and the test requires to be small. If it is less than 0.2 or 0.3, then stop sifting the procedure, and set and is recognized as the first IMF.

After the sifting processing which demonstrates in Figure 4, we can finally get an amount of decomposition IMF component and a residual:The original data can finally be expressed as a sum of IMFs and a residual:where denotes the number of IMFs, denotes the IMF, and is the final residual.

However, EMD suffers from the mode mixing problem which makes EMD fail to represent the characteristics of the original data. Therefore, to overcome the shortcoming of EMD, an additional step adding white noise into the original data is adopted in EEMD model to help extract true IMFs in the data [46].

The procedure of EEMD is presented as follows.

*Step (a)*. Adding a white noise to the original data .

*Step (b)*. Decomposing the data plus the white noise into IMFs through EMD algorithm.

*Step (c)*. Repeating the above two steps, but each time adding different scales of white noise, and then calculating the ensemble means of corresponding IMFs as the final IMFs.

Wu and Huang [46] proved a well-established statistical rule in (16) to control the effect of the added noise:where is the number of ensemble members, is the amplitude of the added white noise, and is the standard deviation of error, which is defined as the difference between the input signal and the corresponding IMFs.

It is suggested that an ensemble of a few hundred and the amplitude of the added white noise fix at 0.2 times the standard deviation of the original signal will lead to an exact result. Unfortunately, there is no literature explicating how to set the best amplitude of the added white noise until now. We set the number of ensemble members as 100 and select the optimal standard deviation of white noise series from 0.1 to 0.2 using -fold cross-validation method. Figure 5 shows the original observed wind speed time series and that by EEMD. And in this study, we remove the first two IMFs with high frequency.

##### 3.3. Cuckoo Search (CS)

Cuckoo search (CS) is a new metaheuristic optimization method firstly proposed by Yang and Deb in 2009 [47]. CS is a new optimization algorithm with an evolutionary process, the basic idea of which is based on the obligate brood parasitic behavior of cuckoos and their way in egg laying and breeding. Actually, the group behavior of searching food is effectively recognized as a random walk because the next move is based on the current location and the transition probability to the next location.

In the survive competition, cuckoos make two main operations: the first is a random search based on the probability of being discovered by a host bird; and the second is a direct search based on Levy flights. With the combination of two operations, CS shows more effective than other optimization algorithms for nonconvex and complex optimization problems.

In the CS method, each nest represents a candidate solution and the best nest is the best solution. For a dimensional optimization problem, each nest is an array with parameters, which represents the current living position of the cuckoo. The main steps of CS are described as follows.

*Step (a)*. Initializing the population : suppose the number of host nests is and the optimization problem has parameters (the parameters in T-S fuzzy networks can be viewed as a nest of optimization search). Then a population of nests is denoted by

*Step (b)*. Updating the new solution through Levy flights: the best path of the Levy flights is estimated by Mantegna’s algorithm. The new solution of each nest is updated by the formulations:where is the individual best, is the global best of the population, is the iteration step size, and rand_{1} is a random number uniformly distributed in the interval from 0 to 1; is defined as , and and have normal distribution with standard deviation and given as follows:where is the distribution factor between 0.3 and 1.99 and is the gamma function.

*Step (c)*. Being discovered by the host bird and randomization of nests: suppose that the probability a cuckoo egg to be discovered by a host bird in its nest is . The action of this discovery also creates a new solution which is defined as follows:where rand_{2} is a uniformly distributed random number and ranges from 0 to 1. and are the random perturbation of the current best nest , respectively. is a number determined based on :After the three steps are done, it is time to stop the algorithm when the maximum number of iterations is reached.

##### 3.4. The Hybrid Model

The main hybrid CS-EEMD-FNN model is proposed to forecast the short-term wind speed in Shandong wind farm in China, which is firstly employing the EEMD method to remove the noise embedded in the wind speed series and then combining FNN with CS. The second hybrid short-term wind speed forecasting model combines FNN with the EEMD approach, that is, EEMD-FNN model. The third hybrid model is based on FNN and the optimized CS algorithm, that is, CS-FNN model. The last one is a single FNN model which is used to compare the performance of the hybrid models. In this paper, the following four forecasting approaches are analyzed:(a)FNN model;(b)EEMD-FNN model;(c)CS-FNN model;(d)CS-EEMD-FNN model.

#### 4. Statistics Measures to Determine the Accuracy of the Forecast

In order to evaluate the performance of the proposed models, comparisons of these forecasted outcomes with the observed wind speed are processed using some statistical evaluation indices. They are the mean absolute error (MAE), squared mean square error (RMSE), and mean absolute percentage error (MAPE) which are defined aswhere is the number of periods of time, and denotes the forecasted wind speed for the time period , and is the actual observation at the same time-point, respectively.

Then these commonly used performance metrics are defined as follows: MAE and RMSE are used to evaluate the absolute error between the observed and the forecasted concentrations. Including the power term in the RMSE makes it more sensitive to excessive values than MAE. In addition to the evaluation of the absolute error, the relative error is estimated using MAPE. For an absolute perfect model, the values of MAE, RMSE, and MAPE are equal to 0.

#### 5. Forecasting Results and Comparative Analysis

To estimate the performance of the proposed hybrid CS-EEMD-FNN forecasting model, three other models, which are the random combination of any two models, are adopted to make one-step ahead predictions. Then in order to evaluate the hybrid CS-EEMD-FNN forecasting model more intensively, multistep ahead predictions are done.

This section shows the forecasting results obtained from FNN model, EEMD-FNN model, CS-FNN model, and CS-EEMD-FNN model for one-step (10 minutes), two-step (20 minutes), three-step (30 minutes), four-step (40 minutes), five-step (50 minutes), and six-step (1 hour) ahead wind speed forecasts.

Once the parameters of FNN are determined, the forecasting results can be obtained by the FNN model with the optimized parameters. A randomly chosen FNN model denoted by FNN model is developed to compare with the optimized FNN models based on the three criteria: MAE, MAPE, and RMSE, shown in Table 1. The forecasting results of one-step, two-step, three-step, four-step, five-step, and six-step ahead forecasting errors of January are shown in Table 2.

Figures 6-7 and Tables 1–3 reveal the following.(i)The more the step of forecast, the lower the accuracy.(ii)These forecasting models adopted in this study are proved to be effective to forecast the short-term wind speed; all the forecasting models get satisfactory performance. Taking the one-step ahead prediction, for example, the MAPE errors of the general FNN, EEMD-FNN, CS-FNN, and CS-EEMD-FNN are 9.12%, 6.46%, 6.47%, and 5.02%, in Jan., respectively.(iii)When comparing the CS-EEMD-FNN model with the EEMD-FNN model, the CS improved the accuracy of multistep ahead forecast of the EEMD-FNN model, and the promoted percentage of MAPE is 4.31%, 5.28%, 5.9%, 6.77%, 7.63%, and 8.64% from one-step ahead to six-step ahead prediction of WTG 25 in Jan., respectively. The CS optimized the initial weights and subsequently improved the global search ability of the EEMD-FNN model. This can accelerate the learning and training process and avoid converging on a local optimal value.(iv)When comparing the CS-FNN model with the FNN model, the CS improved the accuracy of multistep ahead forecast of the FNN model, and the promoted percentage of MAPE is 5.46%, 6.8%, 7.41%, 8.19%, 8.79%, and 9.56% from one-step ahead to six-step ahead prediction of WTG 25 in Jan., respectively.(v)When comparing the EEMD-FNN model with the FNN model, the EEMD decomposition improved the accuracy of the multistep ahead forecast of the FNN model, and the promoted percentage of MAPE is 5.74%, 6.67%, 7.27%, 8.08%, 8.89%, and 9.61% from one-step ahead to six-step ahead prediction of WTG 25 in Jan., respectively. Owing to the EEMD decomposition, the EEMD-FNN model has an excellent ability in signal tracking, which can deal with the wind speed jumping data better.(vi)The EEMD decomposition part in the optimized CS-EEMD-FNN model makes more contribution than the CS part.

#### 6. Discussion and Conclusions

Short-term wind speed forecasting is of crucial importance to improve the stability of grid-connected wind energy system and to avoid unfavorable impacts on the electric network. An accurate short-term wind speed forecasting is thus prerequisite for planning and operating wind power system. This paper proposed a hybrid CS-EEMD-FNN model to the short-term wind speed prediction.

Because of the uncertainty and intermittency, wind speed contains a mass of noisy information, which makes wind speed extremely hard to forecast. As a result, the redundancy information must be removed for good performance of the forecasting model. EEMD is thus used to remove the noise embedded in the wind speed series. In the FNN model, the outputs at any moment only depend on the parameters at that moment. Selecting the optimal parameters is another way to improve the predictive ability of FNN model. With fast convergence rate and low computational complexity, CS is used to tune the optimal parameters in FNN model.

The hybrid CS-EEMD-FNN method is performed on wind speed series in Shandong wind farm, China. The forecast horizons are from 10 minutes to 1 hour. The hybrid CS-EEMD-FNN model is compared to single FNN model, EEMD-FNN model, and CS-FNN model.

Among the hybrid methods proposed in this paper, the CS-EEMD-FNN method has the best performance. Removing redundancy information in the wind speed series is the main reason of its good performance.

Then, the EEMD and FNN parts use the efficient information to build the forecasting model. Although the performance of EEMD-FNN model is slightly inferior to that of CS-EEMD-FNN, it is also satisfactory considering a lower forecasting accuracy. Besides the CS-EEMD-FNN and EEMD-FNN models, the hybrid CS-FNN model also gives satisfactory forecasting outcomes.

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgments

The work was supported by the National Natural Science Foundation of China (Grant no. 71171102).