Abstract
The diagnosis procedure is performed by integrating three steps: multidomain modeling, event identification, and failure event classification. Multidomain model can describe the normal and fault behaviors of hybrid systems efficiently and can meet the diagnosis requirements of hybrid systems. Then the multidomain model is used to simulate and obtain responses under different failure events; the responses are further utilized as a priori information when training the event identification library. Finally, a brushless DC motor is selected as the study case. The experimental result indicates that the proposed method could identify the known and unknown failure events of the studied system. In particular, for a system with less response information under a failure event, the accuracy of diagnosis seems to be higher. The presented method integrates the advantages of current quantitative and qualitative diagnostic procedures and can distinguish between failures caused by parametric and abrupt structure faults. Another advantage of our method is that it can remember unknown failure types and automatically extend the adaptive resonance theory neural network library, which is extremely useful for complex hybrid systems.
1. Introduction
Hybrid systems include both continuous variable dynamic systems (CVDSs) and discrete event dynamic systems (DEDS) [1] and consist of components with continuous and discrete behavior. Because industrial systems have become more complex, most have the continuous and discrete features of hybrid systems and should be described simultaneously by differential and difference equations. There has been much research into the modeling, analysis, and design of hybrid systems, and there have been many achievements in this field [2].
Generally, procedures for diagnosing hybrid systems are designed using extended observers that estimate the system state under normal and fault conditions [3–5]. Researchers describe systems using hybrid automata [6–8], Petri nets [9–12], or switched systems [13, 14] and make a diagnosis according to values measured from a real system. Alternatively, they design an observer [15–19] to describe a system [19], estimate a system’s status (failure or not) [18], or improve the stability [14]. Fault detection is useful and necessary for normal linear and nonlinear system when designing system controller with high robust ability [20, 21]. However, it is difficult to develop a procedure for diagnosing hybrid systems with complex continuous and discrete dynamics, because quantizing the complex behaviors using a finite set of states and events results in large nondeterministic models, degrades the performance, and increases the computational requirements of the diagnostic algorithms.
Daigle et al. [22] presented a novel approach to construct diagnostic procedures for isolating single, abrupt faults in hybrid systems. Their technique was based on a qualitative abstraction of measured deviations from normal behavior, whereas the procedure is not applicable to the hybrid system with several different CVDS. A fault classification framework for a hybrid system has been created by Tromp et al. [9], the hybrid system was modeled based on a Petri net, and the system normal behavior and fault behavior were described as certain discrete event and random discrete event separately. The author further performed the fault diagnosis by obtaining a system event transition sequence based on the measured values from the real system. Philippot et al. [15] developed a discrete event model-based approach for fault detection and isolation of manufacturing systems. Alexandre considers a system as a set of independent plant elements (PEs), each being composed of a set of interrelated parts of plant (PoPs) modeled by a Moore automaton. The PoP model was affected only by its local behavior. Degraded and faulty behaviors were added to each PoP model in order to obtain extended behavior.
Zhao et al. [14] established a diagnostic system architecture that integrated timed Petri net modeling and diagnostic components. The author stated that the advantages of model-driven signature analysis are that it can significantly reduce the computational cost of mode estimation and enabled the diagnostic system to operate online. Bhowal et al. [23] proposed an automaton modeling method for discrete time hybrid system. They utilized a variable-based process model where each activity state represented a particular dynamic of the system. A diagnostic construction method was then formulated based on both the measurable and unmeasurable partitions of the variables and relationships with event transitions.
When designing a model-based diagnostic procedure for a hybrid system, it is important to integrate failure behavior modeling, the acquisition and processing of sensor signals, and fault mode. Researchers need to address diagnostic problems at the interface between model-based diagnostic techniques and signature analysis, so that incipient and abrupt faults can be efficiently detected and isolated. There are, however, few specific reports in the literature distinguishing between the different CVDSs of hybrid systems. Most of the current diagnostic algorithms consider variables which have little impact on the hybrid system; so the algorithms are too complex and that their applications are limited [24, 25].
Diagnosis-oriented modeling is essential for a hybrid system, because it is essential for describing all the important factors, such as continuous variable dynamic processes, discrete event dynamic processes, fault behaviors, and uncertain parameters. The rules of hybrid automata modeling are simple but result in an oversized model when describing a complex system. Petri net model design requires an in-depth understanding of the fault mechanisms of components and their transmission characteristics and requires artificially set failure probabilities. Particle filtering will consume considerable computing resources when making diagnostic procedures for a complex system. As with Petri nets, the transition probabilities in Markov chain models should be set manually.
Modelica is one of the most advanced languages for hybrid system modeling. The equation-based object-oriented methodology supported by Modelica provides several advantages when compared with the block diagram modeling approach. Noncausal description is one advantage of models that use differential, algebraic, and discrete equations. Modelica is an object-oriented multidomain modeling language, which can effectively describe hybrid systems [26–30], and CVDS and DEDS simultaneously. Sanz et al. [26, 29] used a free Modelica library (DEVSLib) that supports the parallel DEVS formalism and is mainly designed to model discrete event systems. They developed an object-oriented modeling methodology, supported by Modelica, for multidomain and multiformalism hybrid control systems. Modelica has also been used successfully in multifield modeling for dynamic systems [31, 32].
The adaptive resonance theory neural network 2 (ART2) [33, 34] has various network structures and strong self-learning ability. ART2 can memorize unknown modes through nonsupervised learning ability. In addition, ART2 can be used to create a self-learning function for a diagnostic system when applied to a poor fault a priori information system.
The Brushless DC motor (BLDCM) has the advantages of high efficiency, high starting torque, and high power density and has been widely used in industrial systems with high performance requirements. The BLDCM switches its three-phase armatures on-off based on position feedback signals. The switching changes the topological structure of the power unit of the BLDCM. Besides this, some continuous variables, such as the current of the armature, torque, and rotation speed, vary in real time. Therefore, it is clear that the BLDCM has both continuous variables (such as phase current, output torque, and output rotation speed) and discrete variables (such as the on-off state of the power unit), and it is a typical hybrid system. Moseler and Isermann [35] proposed a model-based method for the detection of faults in BLDCM. The model was derived from the bridge supply voltage/current and the rotor velocity, and the operating procedure was simplified as a single continuous process, whereas the diagnostic procedure needs to be executed with the parameter estimation results and did not take the six operation stages difference into account. Yu et al. [36] presented a diagnostic approach for a BLDCM based on wavelet transforms and artificial neural networks. The difference with Moseler’s method is that the model considers the effects of six different operating sequences through three fault models, the stator winding interturn short circuit fault model, open-switch fault model, and open-winding fault model.
The aim of the present paper is to provide a diagnostic method for hybrid systems. First, multidomain modeling methods are studied to meet diagnostic requirements of hybrid system. Furthermore, the simulation is executed to obtain the response under failure, and the response would be a priori information for training the ART2 library. A hierarchical fault diagnostic engine is then developed based on the characteristics of the hybrid system. Using the event identification results and estimated parameters, the diagnostic engine is designed to achieve the event transition sequence of the dynamic processes of the system. Both the event transition sequence and system estimated parameters are foundation of primary and detailed classification. Primary classification is performed to distinguish the failure classes of the hybrid system, and the detailed classification is used to achieve the exact fault location; an ART2 structure is designed throughout both primary and detailed classification.
2. Fault Diagnostic Methods for a Hybrid System
A hybrid system fault diagnostic scheme is shown in Figure 1. The scheme consists of three parts: multidomain modeling of the hybrid system, fault diagnostic system training (neural network training), and a fault diagnostic procedure.
Multidomain modeling is the first step and needs an in-depth knowledge of the topological structure of the hybrid system. A component library is necessary for integrating multidisciplinary components in a complex system. Fortunately, abundant freely available libraries for Modelica have been developed by early researchers, which make it easier to model hybrid systems. Modeling, simulation, and analysis are fundamental to the diagnosis. Simulation results are used to train the neural network on the diagnostic procedure.
The neural network is trained as follows. Firstly, the system topology and information regarding the system components are used to create a multidomain model of the hybrid. After establishing a clear understanding of all normal and failure behaviors of the system, simulations are used to determine the systems responses under normal and all known fault behaviors. The simulation results are a priori information for subsequent fault diagnosis. The event transition sequences of the system’s dynamic processes are determined using event identification, and some of the parameters of the system’s characteristics are obtained using parameter estimation. The ART2 structure for primary classification is trained according to the results of the transition sequences and parameter estimations. All the failure events are primarily classified using ART2, and every primary classification includes one or more failure events. Finally, we design the ART2 for each class. During the detailed classification process, some further parameters (extracted using wavelet transforms) are entered into the ART2 structure. When integrated with the system failure classification information, these parameters are used to train the neural network to memorize the corresponding characteristics of each event of a failure class. The result of this training is an ART2 library for detailed classification.
The fault diagnostic processes are as follows. Firstly, based on the system running results, the event transition sequence is obtained using the event identification method, and several characteristic parameters are obtained through parameter estimation. The event transition sequence and characteristic parameters are entered into ART2 for primary classification, and the failure classification number is then known. Once we have extracted the ART2 for this class, the system event transition and further characteristic parameters are entered into the detailed classification network. Classification results should be identified using the current network library to distinguish whether the failure event has been previously memorized. A new failure event is then saved into the current ART2 structure for future detailed classification.
2.1. Multidomain Modeling
Hybrid systems include characteristics of CVDS and DEDS; so it is important that model-based fault diagnostic technique can exactly describe hybrid systems. A hybrid system model must be established and then applied into an actual diagnostic procedure. Modeling should consider a system’s normal behavior, fault behavior, random fluctuation of parameters, and environmental perturbations.
The behavioral modes of the components in a hybrid system can be expressed as follows.
Assuming the system has components, the th component has several behavioral modes and the th behavioral mode is described as the following dynamic system with input and output variables: where is the state variable vector, is the input variable vector, is the output variable vector, and , , and are the structure parameters of the th behavior of the th component.
Behavior modeling based on Modelica is an object-oriented modeling method. First, the hybrid system is divided into several subsystems domain, and interfaces are created between different domains. The behavioral model of the component in each subsystem is then expressed based using (1). The behavior of each component includes normal and fault behaviors. The model structure is shown in Figure 2.
2.2. Theories and Methods for Neural Network Training
This section introduces some theories and methods related to fault diagnostic system training (neural network training), including the event identification method, parameter estimation method, and ART2 structure.
2.2.1. Event Identification
Assuming the hybrid system has normal events, represents the continuous variable dynamic process of the event . The transition path is where represents the transition process from event to event and , is the initial event.
According to (2), when , the system will be within the CVDS.
All the failure events can be separated into two categories. One is a transition failure: the ideal CVDS after transition time should be the process , and a transition failure has occurred when the real CVDS is not the same as process . The other category is a nontransition failure: the system does not transit as it should do between and (), which causes the occurrence of some other process that the system should not perform.
The failure class is primarily based on the event transition sequence of the hybrid system. The real event transition sequence of the system is obtained using event identification; then, the real event transition sequence is compared to the ideal one. Failure occurs if they are not identical. Moreover, primary classification can distinguish different failure classes and positions using event transition sequences.
2.2.2. Parameter Estimation Method
The parameter estimation method can be used to evaluate system parameter values for the mathematical model using an optimization method and the input and output data of the system. The results are useful when diagnosing the system and are also necessary when training the neural network for primary classification.
The hybrid system of (1) with a single output can be described by a differential equation: where is the input vector, is the state vector, and and are the system structure parameters (which can be evaluated using the least mean square (LMS) error algorithm).
The single output in (3) can be solved analytically in a linear time invariant system: where is a time variant function relevant to the real-time input and state vector and is the structure parameters vector
The equation error has been introduced for estimating parameters [35, 37]
Assuming the sample time is , is the number of sampling data sets, and , and the LMS algorithm uses to express the equation error in discrete time
Therefore, the structure parameters can be estimated using where is the step size.
A hybrid system has many failures with different characteristics; so failures cannot be diagnosed by only estimating the structural parameters of the system model using parameter estimation. More powerful signal processing technologies should be used to process the output signal and are necessary for training the neural network structure using characteristic parameters of the hybrid system for detailed classification.
The wavelet transform method is used to extract the hidden characteristic parameters of a hybrid system. For any function of time , the continuous wavelet transform is defined as where , , , is a scale parameter, is a displacement parameter, is the wavelet transform of the function , is the wavelet base function, and is the complex conjugate function of .
For the wavelet transform, the most important process is to find an appropriate base function according to the characteristics of the signal and other specific requirements.
In this research, the output data of the hybrid system is processed using wavelet transforms. By changing the scale parameter , we can discover features of the signal. Thus, the characteristic parameters of the output signal can be extracted and used as input to the neural network training for detailed classification.
2.2.3. ART2
The overall structure of the ART2 was created using adaptive resonance theory [33]. Long term memory from down to up (LTM_U) was integrated into the response layer. The filter and orienting neurons were added to identify failures more easily. The ART2 structure is shown in Figure 3.
The ART2 structure includes six parts: sensor layer, response layer, filter neuron, orienting neuron, LTM_D, and output layer. The weighting vectors in the response layer ( are the corresponding weights of the neuron in the same response layer) and the output vector in LTM_D are the center reference vectors of the from down to up and from up to down memory modes. An angle threshold is set as the classification criterion. By calculating the angle between the two vectors of a new mode and memorizing the mode in the network, it is possible to judge whether the two modes are in the same class. The angle threshold is expressed using a vigilance parameter, which is saved in a neuron. The problem caused by uncertain a priori information can be solved by adjusting the vigilance parameter. ART2 is able to distinguish a new mode and memorize it; thus a memory function is designed to memorize and save the new failure mode. When diagnosing the hybrid system, the mode can be loaded; therefore, the neural network is continually improved and its diagnostic abilities improve over time.
2.3. Failure Event Classification Based on the Neural Network
In the framework of presenting hierarchical fault diagnosis engine, event identification results tend to be an important part for obtaining the event transition sequence of hybrid system dynamic processes. Failure event classification is divided to two separate processes, primary and detailed classification. Primary classification is designed to distinguish the failure classes of the hybrid system. Detailed classification is implemented using an ART2 artificial neural network and is able to precisely identify the fault location.
2.3.1. Primary Classification
We assume that the hybrid system has types of failure modes and define the failure mode set
We complete the primary classification by integrating the event transition sequence using event identification and the parameter estimation results. Assuming that all failure modes are divided into subclasses, we generate the primary classification rules based on the event transition sequence which can be used to diagnose a real system where and .
To simplify the diagnostic method, every failure mode subclass, , should be marked accordingly. The corresponding characteristic parameters for every subclass are extracted during the diagnostic process.
We designed an ART2 structure with inputs and 1 output for the primary classification, as shown in Figure 4. Here the inputs are the characteristic parameters and (or) the event transition sequences of the hybrid system. A simple perceptron is connected to the ART2 output, which translates the output to the failure class number.
2.3.2. Detailed Classification
After the primary classification, each failure class must be diagnosed in detail and memorized by ART2. Each failure subclass has its own ART2 structure for detailed classification. The characteristic parameters of each class are used to train its distinct ART2 structure. The trained neural networks should be able to distinguish different failure events of each class.
We have designed an ART2 structure for detailed classification according to the requirements of hybrid system diagnostics, as shown in Figure 5. The inputs are the characteristic parameters of some failure event, which may be the result of parameter estimation or wavelet transformation (see Section 2.2.1). A simple perceptron is connected to the ART2 output, which translates the output to a failure event number.
3. Study Case: BLDCM Fault Diagnosis
3.1. BLDCM Introduction
We select a three-phase wye-connected BLDCM as the study case. The structure of the three-phase full bridge-driven circuit of the BLDCM is shown in Figure 6. VT1, VT2, VT3, VT4, VT5, and VT6 are power units (MOSFET) switched by the output signal of a HALL sensor.
At any moment, only two power units are in the ON state simultaneously: one is in the upper bridge and the other in the lower bridge. Each power switch works continuously within 120 electrical degrees, and the phase on-off state changes every 60 electrical degrees. The on-off shifting rule is shown in Figure 7.
Through analyzing the BLDCM structure, it is obvious that the BLDCM power units have six logic states. Their corresponding equivalent circuits are as shown in Figure 8.
Figure 8 shows that the pairs of states 1 and 4, 2 and 5, and 3 and 6 are similar, except for opposite currents in the armatures. The equivalent mathematical models of State 1 , State 2 , and State 3 are shown in the following equation: where is the armature voltage of the th state, the bus current of the th state, the motor rotation speed of the th () state, the armature resistance of phase , the armature resistance of phase , the armature resistance of phase , the armature inductance of phase , the armature inductance of phase , the armature inductance of phase , and the EMF constant of the th state.
3.2. BLDCM Multidomain Modeling and Simulation
We modeled the BLDCM and its drive circuit model using the Modelica language, as shown on the right in Figure 9. It consists of a three-phase armature, back EMF module, position sensor, and inertia module. The two parts of the system are integrated by an interface between their domains.
The parameters of the BLDCM are bus voltage , rated speed r/min, stator equivalent resistance Ω, stator equivalent inductance H, back EMF constant Nm/A, pole pairs , rotor inertia moment kgm2, and load Nm.
The modeling and simulating results for normal behaviors and some common failure events are shown in Figure 10.
(a) Rotation speed
(b) Bus current
(c) Armature current
Figure 10(c) shows that when VT5 breaks, the armature can only be ON in a negative direction; so the armature has only negative current. One power unit fault will affect two of the six states during one cycle. The bus current of these two states is nearly zero (Figure 10(b)), which causes the motor rotation speed to decrease and oscillate simultaneously.
As shown in Figure 10(c), when the armature breaks, it has zero current, and it is obvious it will affect four of the six states during one cycle. The bus current of these four states is nearly zero. This causes the motor rotation speed to decrease with strong oscillations.
3.3. BLDCM Fault Diagnosis
3.3.1. Failure Event Analysis
BLDCM failure events can be categorized into four main classes: motor power unit (MOSFET) break, armature break, interarmatures shorting, and armature resistance increase. This paper focuses on these four classes of failure events and diagnoses. We consider the eighteen failure events shown in Table 1.
3.3.2. BLDCM Event Analysis and Identification
As described above (Figure 7), only two power units work simultaneously, and the phase condition changes every one-sixth cycle (60 electrical degrees). Six events occur in sequence during one cycle, and the ON sequence of the armatures is A-B, A-C, B-C, B-A, C-A, and C-B. The on-off control of the MOSFETs and armatures is related to the rotator position, as shown in Table 2 (MOSFET: 1 ON, 0 OFF; armature: 1 positive conduction, −1 positive conduction, 0 OFF).
The analysis in Table 2 shows that the BLDCM has six main events, and the event transition sequence is as shown in Figure 11.
One cycle completes within 360 electrical degrees. According to Table 2, a BLDCM system event can be identified by analyzing the control signal of six MOSFETs. Furthermore, one cycle of bus current information can be extracted and be used for further diagnosis by parameter estimation and data preprocessing.
3.3.3. Parameter Estimation
The following variables are used for parameter estimation:
Combining (12) with (13) results in the following:
Equations (14) can be generalized as
Based on the requirements of the parameter estimation method in Section 2.2.1, (15) can be rewritten as
After sampling the data of the bus voltage , bus current , and rotation speed , estimating the values of , , , , , , and using the relationships expressed in (13), we obtain where , , and are the estimated values of the BLDCM resistances of the three armatures, , , and are the estimated values of the inductances (self-inductance + mutual inductance) of the armatures, and is the estimated value of the back EMF constant. Defining and , the BLDCM estimated values under different failure events are shown in Table 3.
As shown in Figure 10, both MOSFET and armature failures will cause the bus current to go to zero. Therefore, only two variables are defined: one is the current threshold and the percentage % which represents the percentage of bus current exceeding during one cycle. Using simulations and analyses of all the failure events, we calculate the % as shown in Table 4.
3.3.4. Training the Neural Network Structure for Primary Classification
Based on the value of %, , and , all 18 failure events can be divided into five primary classes: normal, armature resistance increase, MOSFET break, armature break, and interarmatures short. We designed an ART2 structure with three inputs and one output for primary classification, as shown in Figure 12. A simple perceptron is connected to the ART2 output, which translates the output to the failure class number. As the inputs to the ART2 structure, parameters %, , and are used to train the neural network automatic classifications. In Figure 12,
The output of the simple perceptron is the failure class number. The input/output relationships of the neural network designed for BLDCM primary classification are shown in Table 5.
3.3.5. Training the Neural Network Structure for Detailed Classification
The primary classification only determines the information of the failure class; so exact information, such as failure event position and mechanism, must be obtained by further analysis of the bus current. In this research, we have selected one cycle of data for the bus current to be processed using the wavelet transform (9). We use this to find the varying adjustments of the wavelet coefficient under different scale parameters, which are then used for detailed classification of the BLDCM.
After the primary classification, three failure classes of P2, P3, and P6 can be determined by using the fluctuations of the estimated value of the resistance of all three armatures. The two failure classes of P4 and P5 can be determined using the wavelet coefficients of the six states during one cycle. Owing to the characteristics of the BLDCM bus current, we determine that the appropriate wavelet base function was the Gauss wavelet (shown in (19))
For the three failure classes of P2, P3, and P6, we designed three ART2 structures with three inputs and one output (Figure 5). The three inputs are , , and
A priori information (input and output) of all three ART2 structures is shown in Table 6.
For the two failure classes of P4 and P5, it is difficult to perform a detailed classification using only the estimated resistance and inductance values. The wavelet coefficients of all six states were obtained by processing the bus current data using the continuous wavelet transform method. Bus currents under failures P4 and P5 are shown in Figure 13.
The data are processed using the Gauss2 wavelet (see (9)). Figure 14 shows the wavelet coefficient curves obtained after extracting the wavelet coefficient (scale parameter ). The mean value of the wavelet coefficients is extracted as the characteristic parameters (see Table 7), which are then used for further diagnosis.
We designed two ART2 structures with six inputs and one output for the two failure classes of P4 and P5, as shown in Figure 15. The inputs are the wavelet coefficient mean values of the six states (). These six ART2 structures were trained using the wavelet coefficients of the six states for detailed classification.
A priori information for training the P4 and P5 class failures is shown in Table 7.
3.3.6. Validation of the Fault Diagnostic System
After training the neural network, the diagnostic engine has the ability to diagnose the BLDCM. We introduced a fault into the BLDCM system: a series connecting resistance ( Ω) with the armature B. The measured rotation speed and bus current are shown in Figure 16.
(a) Rotation speed
(b) Bus current
By measuring the position signal using a HALL sensor, the event identification can be finished to obtain the event transition sequence. Once the system failure occurs, three cycles of bus current are collected and used for parameter estimation. The results of the parameter estimation are shown in Table 8.
Based on the bus current shown in Figure 16, % is 100%. Combining this with the estimated parameters of and , the inputs and output of the ART2 for primary classification are shown in Table 8, the primary classification output number is P2, and the failure is in class B1.
We obtained the detailed classification parameters for class B1 (see Table 6) using parameter estimation (Table 8). These were entered into the detailed ART2 and obtained an output number of 2-2. According to Table 6, the failure event number is , which means the phase B armature resistance has increased (see Table 1). This is consistent with the actual failure injected into the system.
4. Conclusions
This research developed a hierarchical fault diagnostic system, which integrated the theories on multidomain modeling, event identification, and ART2s. We demonstrated the validity of the fault diagnostic system using a typical hybrid system: BLDCM.
The three main achievements of this research are as follows.
The normal and fault behaviors of the hybrid system appear to be able to be modeled and simulated conveniently by the Modelica language, which would provide sufficient prior information for the fault diagnostic system.
Using event identification, parameter estimation, wavelet transform, and ART2 technologies, we implemented a hierarchical fault diagnostic engine. Firstly, we obtained the event transition sequence of the system dynamic process using event identification and the estimated parameters of the hybrid system based on the extracted data of some subprocesses. The failure event was then primarily classified by integrating the event transition sequences and estimated parameters. Finally, we designed an ART2 structure for each failure class which tends to efficiently finish the detailed classification of the failure event.
We have demonstrated that our fault diagnostic method was able to identify both known and unknown failure events using a priori information of the hybrid system. Particularly, for a system with less response information under a failure event, the accuracy of diagnosis seems to be higher. Moreover, the presented method integrates the advantages of current quantitative and qualitative diagnostic procedures and can distinguish between failures caused by parametric and abrupt structure faults. Another advantage of our method is that it can remember unknown failure types and automatically extend the ART2 library, which is extremely useful for complex hybrid systems.
The ART2 used in this paper memorizes failure information through the connection weights between neurons. Therefore, diagnosis results marked by meaningless numbers have poor interpretability. A proposed future study is to develop a corresponding expert system to create the relationship between the outputs of the ART2 network and real phenomena.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
Research for this paper was partially supported by National Natural Science Foundation of China (Grant no. 51205007) and Ph.D. Foundation of Ministry of Education of China (Grant no. 20121102120046). In addition, the authors would like to thank the anonymous reviewers who have helped to improve the paper.