Abstract
This paper considers a joint optimal design of admission control and resource allocation for multimedia services delivery in highspeed railway (HSR) wireless networks. A stochastic network optimization problem is formulated which aims at maximizing the system utility while stabilizing all transmission queues under the average power constraint. By introducing virtual queues, the original problem is equivalently transformed into a queue stability problem, which can be naturally decomposed into three separate subproblems: utility maximization, admission control, and resource allocation. A thresholdbased admission control strategy is proposed for the admission control subproblem. And a distributed resource allocation scheme is developed for the mixedinteger resource allocation subproblem with guaranteed global optimality. Then a dynamic admission control and resource allocation algorithm is proposed, which is suitable for distributed implementation. Finally, the performance of the proposed algorithm is evaluated by theoretical analysis and numerical simulations under realistic conditions of HSR wireless networks.
1. Introduction
With the rapid development of highspeed railway (HSR) around the world, the wireless communication in HSR networks plays an important role in the recent years [1]. On the one hand, more and more data related with the railway controlling information needs to be transmitted between the train and the ground such that the safety can be guaranteed and the transportation efficiency can be significantly improved. On the other hand, the passengers in the train have an increasingly high demand on multimedia services. However, these requirements on high throughput impose a great challenge over the HSR communication designs due to the fastvarying channel, train penetration loss, and so forth.
There have been some recent works to improve the throughput in HSR wireless networks. A twohop HSR network architecture was proposed in [2] to provide high datarate services. A HSR communication system based on radio over fiber technology was proposed in [3], which can increase the system throughput and help to reduce the number of handoffs. Multiinput multioutput (MIMO) antennas were employed to improve the throughput performance of the HSR wireless networks [4, 5]. However, these works were carried out only to improve the throughput performance in HSR wireless networks. Since the buffering is involved at network devices, for example, content servers, it is necessary to consider not only the throughput performance but also the queue stability in HSR wireless networks.
Admission control and resource allocation, as critical parts of radio resource management, play important roles in improving the throughput and ensuring queue stability. In the literature, the energy constrained control algorithm was proposed in [6] to stabilize the queue and maximize the throughput by Lyapunov optimization theory. Paper [7] studied the joint scheduling and admission control problem in a single user scenario and an online learning algorithm was proposed based on the Markov decision process approach and stochastic control theory. However, these existing schemes designed for general communication systems are not easily extended to the scenario considered in this paper, due to the following reasons: in HSR wireless networks, the channel condition cannot remain at the same level because of the fastvarying distance between the base station and the train, which causes that the power control along the time has a large influence on system transmission performance [8]; many types of services with different qualityofservice (QoS) requirements and priorities should be supported [9], which makes the admission control and resource allocation for multiple services more challenging.
In HSR wireless networks, few studies have been conducted on admission control and resource allocation. A scheduling and resource allocation mechanism was presented in [10] to maximize the service rate in HSR networks with a cell array architecture. In [11], a multidimensional resource allocation strategy was proposed in downlink orthogonal frequencydivision multiplexing (OFDM) system for HSR communications. The optimal resource allocation problem in a cellular/infostation integrated HSR network was investigated in [12], which considered the intermittent network connectivity and multiservice demands. In a relayassisted HSR network, [13] studied delayaware fair downlink service scheduling problem with heterogeneous packet arrivals and delay requirements for the services. Paper [14] proposed an effective admission control scheme considering different service priorities for HSR communications with MIMO antennas. However, to the best of our knowledge, the joint admission control and resource allocation problem under the average power constraint in HSR wireless networks is still an open problem.
The main contribution of this paper is a stochastic optimization framework for transmitting multimedia services in HSR wireless networks, which focuses on the joint admission control and resource allocation problem under the average power constraint. Firstly, the joint admission control and resource allocation problem is formulated as a stochastic optimization problem, and then the problem is transformed into a queue stability problem with the help of virtual queues. By the driftpluspenalty approach [15], the transformed problem can be decomposed into three separate subproblems: utility maximization, admission control, and resource allocation. The former two subproblems are easy to handle and the distributed solutions can be obtained directly, while the mixedinteger resource allocation subproblem is transformed into a single variable problem and a distributed packet loading scheme is developed with guaranteed global optimality. We further propose a dynamic admission control and resource allocation algorithm, which is suitable for distributed implementation in HSR wireless networks. Finally, we present the analysis of algorithm performance by theoretical derivations and simulations under realistic conditions for HSR wireless networks.
1.1. Relation to Prior Work
The Lyapunov drift theory has a long history in the field of discrete stochastic processes and Markov chains [16]. It can be used to directly analyze the characteristics of the control policies in the stochastic stability sense and plays important roles in the dynamic control strategies in queuing networks [17]. Stabilizing queuing networks by minimizing Lyapunov drift was pioneered by Tassiulas and Ephremides in [18]. The Lyapunov drift theory was then extended to the Lyapunov optimization theory [6], which enables optimization of time averages of general network utilities subject to queue stability. A general framework for solving the stochastic network optimization problem based on Lyapunov optimization theory was developed in [15]. This framework has been extended to minimizing a driftpluspenalty expression in [6, 7, 17, 19, 20] for joint queue stability and time average utility optimization. For the engineering applications of Lyapunov optimization theory, interested readers are referred to the aforementioned references for the details.
Our approach in the present paper treats the joint admission control and resource allocation problem associated with average power constraint using Lyapunov drift and Lyapunov optimization theory from [15]. This is the first time, to the best of our knowledge, that the Lyapunov optimization theory is extended into the HSR wireless networks. Considering the features of HSR wireless networks, the Lyapunov optimization theory is successfully applied for solving the joint admission control and resource allocation problem in HSR wireless network.
1.2. Outline of Paper
The rest of the paper is organized as follows. Section 2 describes the system model. The problem formulation and transformation are provided in Section 3. A distributed dynamic admission control and resource allocation algorithm is proposed in Section 4. Some numerical results and discussions are shown in Section 5. Finally, conclusions are drawn in Section 6.
Notations. In this paper, denotes expectation. . and mean the maximum and minimum between and , respectively.
2. System Model
In this paper, a twohop HSR wireless network architecture is considered, as shown in Figure 1, which consists of a backbone network, content servers (CSs), several base stations (BSs), a relay station (RS), and some users in the train. The BSs deployed along the rail line can provide continuous data packets delivery. The distributed CSs connected to the BSs via wireline links are deployed in the backbone network to offload the data traffic [21]. The RS with powerful antennas installed on the top of the train is used for communicating with the BSs so that the large train penetration loss can be well resolved. The RS is further connected to the access points (APs) which can be accessed by the users inside the train. Thus, there is a twohop wireless link, consisting of the BSRS link and the APUsers link. If the users on the train request multimedia services during a trip, the data packets of the requested services are then delivered from the corresponding CS to the RS via a BS. Suppose that the data transmission rate in the APUsers link is sufficiently large; hence the data packet can be successfully received if it has been delivered to the RS.
2.1. TimeDistance Mapping
Consider a train traveling from an origin station to a destination station within the time duration . The whole time is divided into slots of equal duration . Without loss of generality, we assume that the train starts at the centre of the first cell and the train moving speed during the slot keeps constant, denoted by ; thus the traveled distance until slot is given by . The train location between two adjacent BSs at slot is , where is the cell radius. Define a timedistance mapping function as the distance between BS and RS at slot ; that is, , where and is the distance between each BS and the rail line as shown in Figure 1. The mapping function can be expressed by Here we assume that the distance does not change within slot since is very small.
2.2. Physical Layer Model
For HSR wireless networks, the channel condition cannot remain at the same level due to the fastvarying distance between BS and RS. Only the lineofsight (LOS) path in the BSRS link is available at most of the time, which was confirmed by engineering measurements [22, 23]. The service provided by the independent identical distributed (i.i.d.) fading channels is a deterministic timelinear function, just like the AWGN channel [24]. Therefore, the wireless channel in the BSRS link can be assumed to be an additive white Gaussian noise channel (AWGN) with LOS path loss [8]. At the same time, the power control along the travel time in HSR wireless networks is important. Denote by the transmit power of the BS at slot , which is limited by the maximum value and average value . With the help of mapping function and according to Shannon's theorem [25], the transmission rate of the wireless channel between BS and RS at slot can be expressed by where , is the system bandwidth, is the noise power spectral density, and is the path loss exponent. Suppose that the packets have equal size of bits; hence the link capacity at slot can be denoted as the maximum number of packets; that is, . Note that the maximum capacity can be obtained when and .
2.3. Service Model
Assume that there are types of services in the HSR wireless networks and the service type set is denoted by . We further assume that is equipped with a buffer and can provide service , for . Let represent the packet arrival vector, where denotes the number of new arrival packets of service at slot . The packet arrival process for each service is assumed to be i.i.d. across slots. Suppose that, in general, follows a truncated Poisson distribution with average arrival rate for service , and the distribution can be written as where denotes the maximum number of arrival packets per slot for service and can be found assuming .
Let represent the vector of current queue backlogs, where denotes the number of packets at the beginning of slot in the buffer of . The dynamics of each buffer are controlled by admission control (AC) and resource allocation (RA) actions. Specifically, at each slot, the AC action determines the number of packets from the newly arriving packets to be stored into the buffer. And the RA action determines the number of packets removed from the buffer for transmission. Let and denote the AC action and RA action for service at slot , respectively. Thus, the queue dynamics can be characterized by Notice that for any slot , without AC actions, . Here we assume that the arrival packets at slot can only be transmitted at slot .
3. Problem Formulation and Transformation
3.1. Problem Formulation
In this paper, the objective of the joint AC and RA problem is to maximize a sum of utility functions under time average constraints by designing a dynamic algorithm over a trip of the train. We define that is a utility function to present throughput benefit for service , which is nondecreasing concave continuous with . Throughout this work, the following notation for the longterm time average expectation of any quantity is defined: In particular, represents the average queue backlog in the buffer of and represents the average power consumption along the travel time. Here we introduce the definitions of queue stability as follows [15].
Definition 1. A single queue is mean rate stable if .
Definition 2. A single queue is strongly stable if .
From Definition 2, a queue is strongly stable if it has a bounded time average backlog. Strong stability implies mean rate stability according to [15]. Throughout this paper, we use the term “stability” to refer to strong stability. Define as the observed network event at slot . For each slot , observing the event and the queue state , the AC actions and RA actions should be made for . The joint AC and RA problem is formulated aswhere (6b) corresponds to the power constraint and (6c) corresponds to the queue stability constraints for all queues. Problem (P1) is a stochastic optimization problem [15], but it cannot be solved efficiently owing to the difficulty from the objective function (6a) and the average power constraint in (6b). In order to better characterize the problem () and develop an efficient algorithm, we consider the problem transformation, which consists of two steps, that is, objective function transformation and average power constraint transformation as presented in the following subsections.
3.2. Objective Function Transformation
Since problem (P1) involves maximizing a function of time averages, it is hard to handle. Based on the dynamic stochastic optimization theory [15], it can be transformed into an equivalent problem that involves maximizing a single time average of a function. This transformation is achieved through the use of auxiliary variables and corresponding virtual queues with queue evolutions: where the initial condition is assumed that ,. Intuitively, the auxiliary variables can be viewed as the “arrivals” of virtual queues , while can be viewed as the service rate of such virtual queues.
Then, we consider the following transformed problem: Constraint (8b) corresponds to the stability of the virtual queue , since and are regarded as the timeaveraged arrival rate and the timeaveraged service rate for the virtual queue , respectively. Specifically, from (7) we can obtain that . By summing this inequality over time slots and then dividing the result by , we have that . With , taking expectations of both sides yields that . If the virtual queues are mean rate stable, then , so that constraint (8b) can be satisfied. Notice that we will prove the strong stability of the virtual queues in Lemma 7 later.
Lemma 3. Problem (P1) and problem (P2) are equivalent.
Proof. The proof of Lemma 3 follows [26] and a sketch of the proof is provided in Appendix A.
3.3. Average Power Constraint Transformation
To handle the average power constraint in (6b), we define a virtual queue for each , which has the following dynamic update equation: where and can be viewed as the “arrivals” and “offered service” at slot , respectively.
Based on [15, Chapter 4], if the virtual queue is mean rate stable for , that is, , then the average power constraint can be satisfied. This holds because if the backlog in the virtual queue is stabilized, it must be the case that the time average arrival rate (corresponding to ) is not larger than the service rate (corresponding to ). Therefore, the average power constraint in (6b) can be transformed into a single queue stability problem.
4. The Distributed Dynamic AC and RA Algorithm
In this section, the dynamic stochastic optimization approach is applied to solve problem (P2), which seeks to maximize the sum of timeaveraged utility functions subject to queue stability constraints. Firstly, the problem (P2) is decomposed into three separate subproblems by the driftpluspenalty approach. Then a distributed dynamic AC and RA algorithm is proposed. Finally, the performance of the proposed algorithm is analyzed by theoretical derivations.
4.1. Lyapunov Drift
Define and as a vector of all virtual queues and for , respectively. We denote by the combined vector of all virtual queues and all actual queues; namely, The quadratic Lyapunov function is defined as [15] Then the oneslot conditional Lyapunov drift at slot is given by which admits the following lemma.
Lemma 4. Under any AC actions and RA actions at slot , and for any value of , we have where is a finite constant defined by and is defined by
Proof. The proof of Lemma 4 is provided in Appendix B.
4.2. The DriftPlusPenalty Expression
Instead of directly minimizing the upper bound by taking AC actions and RA actions, we desire to jointly stabilize all queues and maximize the sum of utility . The driftpluspenalty theory in [6] approaches this by greedily minimizing the following “driftpluspenalty” expression: where is a parameter that represents the weight on how much we emphasize the sum utility maximization.
We observe that the objective function in (16) is of separable structure, which motivates us to determine the auxiliary variables and AC actions as well as RA actions in an alternative optimization fashion. The overall minimization problem (16) is decomposed into three separate subproblems. Specifically, isolating the variables from (16) gives the following utility maximization subproblem: Similarly, isolating the AC actions from (16) leads to the following admission control subproblem: Also, isolating the RA actions from (16) gives the following resource allocation subproblem: where is related to since a larger requires more power consumption. These separate subproblems can be computed in a decentralized fashion, as stated below.
4.3. Utility Maximization
The utility maximization subproblem ((17a) and (17b)) can be decoupled into separate maximization problems. Specifically, keeps track of and determines the optimum by solving the following problem:
Notice that the key point to solve (20a) and (20b) is the choice of the utility function, which is contingent on the purpose of the networking application or the prerogative of HSR network designer. For example, in order to represent the maximum desired delivery ratio for each service, the piecewise linear utility function can be considered for service as follows: where and represent the priority and the maximum desired delivery ratio of service , respectively. In general, and for . Thus the optimal solution to problem (20a) and (20b) is given by Alternatively, the following strictly concave function can serve as the utility function for service : which can be regarded as an accurate approximation of the proportionally fair utility function if the same is selected with a large value for all . In this case, the optimal solution to problem (20a) and (20b) can be obtained by where the operation is equal to if , if , and if .
4.4. Admission Control
The admission control subproblem ((18a) and (18b)) can be also decoupled into separate maximization problems. Specifically, chooses the AC action by solving the following optimization problem: It is immediate to see that the optimal solution depends on the queue backlog of and , which is given by We note that (26) is a simple thresholdbased admission control strategy. On the one hand, when the queue backlog is not larger than the threshold , then all the newly arriving packets are admitted into the buffer in . Essentially, this not only reduces the value of so as to push closer to , but also increases the average throughput so as to improve the utility. On the other hand, when the queue backlog is larger than the threshold , then all the newly arriving packets will be dropped to ensure the network stability. Finally, we emphasize that the AC actions for all services are made in a distributed manner with only local queue backlog information and packet arrival information.
4.5. Resource Allocation
The resource allocation subproblem ((19a)–(19d)) at slot can be explicitly expressed as The problem (27a)–(27d) is a mixedinteger programming (MIP) problem, including a continuous variable and integer variables , which cannot be solved efficiently [27]. The main difficulty of problem (27a)–(27d) comes from the integer nature of . However, we will show that problem (27a)–(27d) can be transformed into a single variable problem, which is easy to handle. In the sequel of this subsection, we will omit the time index for brevity.
Firstly, as for constraint (27c), when the optimal RA actions are achieved, it can be shown that where . Otherwise we can reduce the value of and such that the objective function can be further maximized without any violation of the constraints in (27b)–(27d). From (28), we have the following power consumption of : and constraints (27b) and (27d) further imply that where .
Secondly, the resource allocation subproblem (27a)–(27d) can be equivalently transformed into a single variable problem as follows: where is given byand is given by with .
Now let us focus on the problem (32a)–(32c) with any given . Clearly, the maximum objective value of (32a)–(32c) can always be achieved by allocating link capacity to the services in the descending order of their backlogs, which is similar to the maxweight algorithm in [15]. Hence, we sort all the services in descending order of and denote the ordered set by . For convenience, we define for , where . One can see that is an increasing function of , and from (30). Therefore the optimal solutions to the problem (32a)–(32c) are given by where such that if ; otherwise for all .
Next, let us focus on problem (31a)(31b). Indeed, we have the following lemma.
Lemma 5. is a unimodal function of over .
Proof. On the one hand, since for , is a monotonically nonincreasing function of . On the other hand, since , is a monotonically increasing function of . Therefore , which is a monotonically decreasing function of . For any , , which implies , so is concave on . Based on [28], if is concave, then is unimodal.
Based on Lemma 5, since is a unimodal function of over , the golden section search method [29] can be used to obtain the global optimal solution to the problem (31a)(31b). However, this method requires the knowledge of all queue backlog information. When the center controller is not available, a distributed resource allocation scheme is highly desirable. Relying on the insights from Lemma 5, we propose a distributed resource allocation scheme, where each CS can communicate with all other CSs, and the network resources are allocated packet by packet.
The proposed packet loading resource allocation scheme is detailed in Algorithm 1. For each slot, each CS exchanges the backlog information with all other CSs and the order of the backlogs is obtained by all CSs (step 2). Then the packets of the services are fetched from the corresponding CS in descending order of their backlogs. When one CS fetches a new packet (step 5), is calculated (step 6). This will be repeated until the optimal condition in step 7 is satisfied, which implies that is the maximum or constraint (30) is violated, and thus step 8 should be performed since the last packet cannot be transmitted. When one CS empties its buffer, it should send the value information of to the next CS (step 12), and the next CS should continue the resource allocation. Here we remark that the optimal solutions can be achieved by the proposed scheme, which can be readily proved by Lemma 5.

4.6. Distributed Dynamic AC and RA Algorithm
Based on the above three separate subproblems, we propose a distributed dynamic AC and RA algorithm as shown in Algorithm 2. All system parameters should be initialized before the trip begins. At each slot, each CS solves three subproblems in steps 4, 5, and 6. At the end of each slot, the queue vector is updated according to (4), (7), and (9). This algorithm will be repeated when the train travels from the origin station to the destination station.
Remark 6 (utilitybacklog tradeoff). Based on [15], the achieved utility differs from optimality by , in the sense that
where is the maximal utility for the problem (P1). It implies that the proposed algorithm can achieve a utility which is arbitrarily close to by increasing . In addition, the actual queue backlog of each service grows linearly with , which is given by
where is a parameter and is defined in (14). Therefore, the above expressions (35) and (36) present a utilitybacklog tradeoff of .
Recalling the utility functions (21) and (23), we observe that they have the maximum right derivatives over the interval . Based on this observation, we obtain the boundedness property on the virtual queue in the following lemma.
Lemma 7. If the utility function has maximum right derivatives , then the backlog of virtual queue satisfies provided that this inequality holds for .
Proof. The proof of Lemma 7 is provided in Appendix C.
5. Numerical Results and Discussions
In this section, we implement the proposed distributed dynamic AC and RA algorithm using MATLAB and present simulation results to illustrate the performance of it. We use the piecewise linear utility functions (21) for all services and summarize the simulation parameters in Table 1. The packet size is set to 240 bits according to [12, 30], and the slot duration is set to 1 ms according to [31]. A single simulation runs the proposed algorithm when the train moves through a cell (30,000 slots).
Figures 2 and 3 explore the throughputbacklog tradeoff with different . In the simulations, we use the same parameters and and different priorities for all services. As shown in Figure 2, the average throughput for each service increases as is increased and the service with high priority gets the large average throughput. Figure 3 presents that the average queue backlogs of all services are linearly increasing with , which demonstrates the behavior in (36). Furthermore, the proposed algorithm can ensure that the average queue backlogs of the services with different priorities are almost the same.
Figure 4 illustrates the achieved delivery ratios for the services with different maximum desired delivery ratios and same arrival rate as well as priority . It can be observed that the large will result in the improvement of the delivery ratio performance. This can be explained as follows: since a larger gives a higher priority on throughput, more packets will be admitted into the buffers, which causes the higher delivery ratio performance. In addition, the delivery ratio for each service is close to its maximum desired delivery ratio when is larger than 40, which implies that the proposed algorithm can archive different maximum desired delivery ratios when a large is chosen.
Figure 5 compares the average power consumption under different arrival rate conditions. In this simulation, we set the same parameters , , and for all services. From this figure, we can see that the average power consumption increases as is increased. This is exactly what happens. A larger results in more packets admitted into buffers, while transmitting these packets will cost more power. As for the same , a larger will cause more power consumption, since more packets will be admitted into buffers based on (26). In addition, the average power consumption can satisfy the average power constraint when the arrival rate is , which is reasonably set in the previous simulations.
Figure 6 describes the backlog update processes of virtual queues for three types of services. In the simulation, we set the same parameters , , and for all services and . From the figure, we can see for all services at all slots, which illustrates the boundedness property on queue backlogs in Lemma 7.
6. Conclusion
In this paper, we formulate the joint admission control and resource allocation problem under average power constraint for multimedia services delivery in HSR wireless networks. With the help of virtual queues, the original stochastic optimization problem is transformed into a queue stability problem, which is decomposed into three separate subproblems by the driftpluspenalty approach. It is worth noting that the optimal solution to the resource allocation subproblem can be obtained by the packet loading resource allocation scheme. Based on the stochastic optimization technique, the dynamic admission control and resource allocation algorithm is proposed, which is suitable for distributed implementation in HSR wireless networks. Furthermore, the performance of the proposed algorithm is analyzed theoretically and validated by numerical simulations under realistic conditions for HSR wireless networks. For future work, we will further investigate the effects of the nonmentioned parts in the communication system, such as frame error check blocks and adaptive channel equalizers.
Appendices
A. Proof of Lemma 3
Let and be the optimal utility of problems (P1) and (P2), respectively.
First, to prove , let be an optimal solution achieving in problem (P2), which includes , , and at slot . Since is concave, based on Jensen's inequality, we have In addition, since the solution satisfies constraint (8b) and is nondecreasing, we further have Since the constraints in problem (P2) include all of the desired constraints of the original problem (P1), is a feasible solution for problem (P1) which gives a utility that is not larger than . Thus we conclude that Next, to prove , let be an optimal solution achieving for problem (P1), which includes and at slot . Since satisfies constraints (6b)–(6f), it also satisfies constraints (8d) of the problem (P2). Further, for all , we set at all time which can satisfy constraints (8b) and (8c). Thus, such choice of together with the solution forms a feasible solution for the problem (P2). By definition, . Therefore, we get From the above analysis, we can conclude that based on (A.3) and (A.4) and that an optimal solution for the problem (P2) can be directly turned into an optimal solution for the problem (P1).
B. Proof of Lemma 4
Recall the evolution equation (9) for the queue and by squaring this equation, we obtainwhere in the final inequality we have used the following facts: and .
Similarly, it can be shown for that Based on (12) and (B.1a)–(B.2), we havewhere is defined by (15) and the last inequality can be obtained as follows. For any slot , any possible packet arrival vector , and any possible as well as RA actions that can be taken, we have where the inequality holds based on , , and as well as (8c), and the equality holds since the constant in the square bracket is independent of queue vector at slot .
C. Proof of Lemma 7
We prove this lemma by induction. Assume that for slot (it holds by assumption at slot ); then we prove it also holds for slot . Firstly, we consider the case . From the queue update equation (7), we can see that this queue can increase by at most at each slot, and thus we have , proving the result for this case.
Secondly, we consider the case . For each slot , decides to maximize the following expression: Based on the property of the maximum derivative, for any , we obtainwhere the equality holds if and only if . Then the algorithm will choose to maximize expression (C.1), and we can obtain Thus, is satisfied for these two cases, which completes the proof.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work is supported by the Fundamental Research Funds for the Central Universities (Grant no. 2014YJS026), the Key Projects of State Key Lab of Rail Traffic Control and Safety (no. RCS2012ZZ004 and no. RCS2010ZT011), the Opening Project of The State Key Laboratory of Integrated Services Networks, Xidian University (Grant no. ISN1409), and the Key Grant Project of Chinese Ministry of Education (no. 313006).