J. Cent. South Univ. (2016) 23: 3248-3255

DOI: 10.1007/s11771-016-3390-8

NSGAII based multi-objective homing trajectory planning of parafoil system

TAO Jin(陶金)^{1, 2}, SUN Qing-lin(孙青林)^{1, 2}, CHEN Zeng-qiang(陈增强)^{1, 2}, HE Ying-ping(贺应平)³

1. Tianjin Key Laboratory of Intelligent Robots, Nankai University, Tianjin 300071, China;

2. College of Computer and Control Engineering, Nankai University, Tianjin 300071, China;

3. Aerospace Life-Support Industries Ltd., Aviation Industry Corporation of China, Xiangyang 441003, China

Central South University Press and Springer-Verlag Berlin Heidelberg 2016

Abstract: Homing trajectory planning is a core task of autonomous homing of parafoil system. This work analyzes and establishes a simplified kinematic mathematical model, and regards the homing trajectory planning problem as a kind of multi-objective optimization problem. Being different from traditional ways of transforming the multi-objective optimization into a single objective optimization by weighting factors, this work applies an improved non-dominated sorting genetic algorithm II (NSGA II) to solve it directly by means of optimizing multi-objective functions simultaneously. In the improved NSGA II, the chaos initialization and a crowding distance based population trimming method were introduced to overcome the prematurity of population, the penalty function was used in handling constraints, and the optimal solution was selected according to the method of fuzzy set theory. Simulation results of three different schemes designed according to various practical engineering requirements show that the improved NSGA II can effectively obtain the Pareto optimal solution set under different weighting with outstanding convergence and stability, and provide a new train of thoughts to design homing trajectory of parafoil system.

Key words: parafoil system; homing trajectory planning; multi-objective optimization; non-dominated sorting genetic algorithm (NSGA); non-uniform b-spline

1 Introduction

A parafoil system is made up of a traditional ram- air parafoil, loads and guidance navigation & control (GNC) equipment. It is a kind of precision aerial delivery system with superior pneumatic performance, excellent gliding ability and easy handling property. It now has been widely used in military, aerospace and civil fields, such as unmanned aerial vehicle (UVA), crew return vehicle (CRV) and guided parafoil air drop system (GPADS) due to its excellent properties [1-3]. With the development of large size ram-air parafoil, global position system (GPS), miniaturized powerful sensor and micro control units (MCUs), the implementation of fully autonomous homing of parafoil system becomes feasible. Homing trajectory planning and optimization is overwhelmingly crucial to realize autonomous precision homing and the pros and cons of homing trajectory affect the effect of homing to a great extent. At present, a variety of homing trajectory planning methods have been recommended, namely the blind-angled planning method, the multi-stage planning method and the optimal planning method, the last of which is now widely concerned by academics [4].

The so-called optimal planning method can be described as searching for a particular trajectory from the initial position to the target that meets specific performance under certain constraints. Such an optimal trajectory planning problem belongs to a category of large-scale complex nonlinear programming, which can be cast as a multi-objective optimization problem. Corresponding solving methods are algorithms based on classical control theory and intelligent algorithms. Scholars both at home and abroad have done a lot of explorations to solve such a problem. Pearson [5] took the lead in carrying out some researches in 1972 and he solved the optimal control problem which could reach the target without any constraints by the differential dynamic programming method. In recent years, especially after 2000, there emerged a large number of research results. Gimadieva [6] used the method of Krylov and Chernous'ko in solving the boundary-value problem for a gliding parachute system. Xiong [4] adopted the conjugate gradient method to solve the problem which is transformed into two-point boundary value problem by the Pontryagin minimum principle. Rademacher et al [7] used optimal control theory to study the onboard trajectory planning of parafoil system. Liu and KONG [8] and Jiao et al [9] respectively adopted the improved particle swarm algorithm in optimizing homing trajectories. Zhang et al [10] and Gao et al [11] used the Gauss pseudo- spectral method in transforming the problem into parameters optimization with a series of algebraic constraints and got a satisfactory solution. Luders [12] presented a new approach based upon chance- constrained rapidly-exploring random trees to online trajectory planning and robust obstacle avoidance for the parafoil terminal guidance problem. Cleminson [13] used a dynamic programming approach in solving the path planning problem of a parafoil. Babu et al [14] adopted the direct multiple shooting method in solving the multi-objective trajectory optimization problem to get the optimal trajectories that minimize performance index functions. Tao et al [15] proposed a homing trajectory planning method of parafoil system based on quantum genetic algorithm.

It can be seen from the aforementioned research achievements that the classical control theory based algorithms including both the indirect and direct methods belong to is traditional search methods based on gradient that are sensitive to initial values. Therefore, the optimization results being good or bad largely depend on speculation of initial values. Intelligent algorithms, such as genetic algorithm (GA) and particle swarm optimization algorithm (PSO), have strong robustness, a wide range of adaptability and are not sensitive to initial values, which make them a better choice for homing trajectory planning of parafoil system. Most of the work has been done trying to solve the multi-objective homing trajectory planning problem by transforming the multi- objective optimization into a single-objective optimization by weighting factors, which is commonly used in multi-objective problem optimization. However, it will bring the following troubles: on one hand, the determination of weighting factors is an optimization process, on the other hand, the final optimization results strongly depend on preferences of weighting factors. The actual operation of choosing weighting factors is usually according to the importance indexes of multi-objective functions, meanwhile the iterative trial and error method is often applied to select a satisfactory solution. In order to overcome the above-mentioned difficulties, the multi- objective optimization algorithms based on Pareto method are introduced to solve the problem. The typical algorithm NSGA II is widely used for its high processing speed and rosy convergence of solutions. This work adopts an improved NSGA II [16-19] to solve the problem of multi-objective homing trajectory planning of parafoil system. Different from existing methods, this algorithm is insensitive to initial values and more able to reflect the nature of multi-objective optimization. The optimization results do not depend on design experiences of weighting factors selecting, and the process avoid the tedious trial and error steps. Evenly distributed multi- objective Pareto optimal solution set under different weighting can be obtained at one operation, which brings more design choices and greater freedom to the practical engineering design.

2 Problem modeling and description

2.1 Mathematical model

The dynamic model of parafoil system is highly complex and nonlinear, which invariably has strong coupling. Therefore, the particle model is quintessentially used in place of the complex high degree of freedom model in homing trajectory planning. By means of this method, calculations process can be greatly simplified.

The wind coordinate system is often adopted in homing trajectory design of parafoil system.The direction of each coordinate axis in the wind coordinate system is consistent with the geodetic coordinate system. The origin of wind coordinate system moves with the airflow and coincides with the origin of geodetic coordinate system at the time when flare landing is implemented. Therefore, the degree and direction of wind can be transformed into the shift position of the origin.

The flight control of parafoil system is realized by manipulations of left or right control ropes on its trailing edge. Jiao et al [20] established a six-degree of freedom dynamic equation of parafoil system. Simulation results of the six degree of freedom model are shown in Tables 1 and 2.

It can be seen from Table 1 that the bilateral lower deviator has little effect on gliding performance of parafoil system except for flare landing. Therefore, theunilateral lower deviator is only considered in control. As shown in Table 2, the increase of unilateral lower deviator leads to the decrease of gliding ratio and roll angle. Consequently, the unilateral lower deviator is usually limited within a certain range for it will damage holistic stability when reducing to a certain degree (when the unilateral lower deviator reaching to 40%, the parafoil system is already in stall condition). Accordingly, in order to simplify the model, the following assumptions are established:

1) The horizontal and vertical velocities remain unchanged under ignoring effects caused by the variation of atmospheric density.

2) The horizontal wind field is fixed and known.

3) The control response of parafoil system is without any delay.

Table 1 Glide ratio with increase of bilateral lower deviator

Table 2 Glide ratio with increase of unilateral lower deviator

Based on the above assumptions and wind coordinate system, the coordinate origin is selected where flare landing is began to implement. The motion equations of parafoil system can be simplified as

                                (1)

where (x, y, z) means the location information of parafoil system; v_s means its horizontal velocity; v_z means its vertical velocity; ψ means its yaw angle; means its yaw rate; u denotes the control quantity, and there exists a one-to-one relationship between u and the unilateral lower deviator.

2.2 Problem description

The problem of homing trajectory planning of parafoil system can be described as the pursuit of control law which can make specific performances optimal under certain dynamic constraints. The requirements and conditions of homing trajectory optimization can be summarized as follows:

1) The landing site is close to the target point.

2) The landing direction is against the wind. This is the necessary condition of flare landing through which the speed of landing can be reduced and the damage of payload can avoid in the process of landing.

3) The feasible manipulation and low energy consumption is required.

According to the particle model and assumptions mentioned above, constraints and objective functions are formulized as follows:

2.2.1 Initial conditions

The initial conditions at the starting time t₀ are known and can be expressed as

                                 (2)

where (x₀, y₀, z₀) denotes the initial coordinate information of parafoil system, and ψ₀ denotes its initial yaw angle.

2.2.2 Terminal constraints

The terminal time of homing is t_f, which can be calculated as t_f =z₀/v_z, and the terminal constraints can be expressed as

                      (3)

where (x_f, y_f, z_f) denotes the coordinate information of the landing site; ζ_wind denotes the direction of horizontal wind. The constraint is to ensure that the landing direction is against the wind.

2.2.3 Control constraint

Due to the limitation of motors, the control quality of motors should be restricted by the following constraint inequality:

                                   (4)

where u_max denotes the allowed maximum control quantity, which also corresponds to the minimum turning radius of parafoil system.

2.2.4 Objective functions

According to the particle model and under the constraints mentioned above, the requirements of homing trajectory of parafoil system are that the landing site is as close as possible to the target; flare landing is necessary for smooth touch down; and the control input is as small as possible to insure the minimal energy consumption. Based on these requirements mentioned above, three objective functions are established as follows:

                    (5)

where J₁ denotes the minimum energy consumption; J₂ denotes the minimum landing error; J₃ denotes the demand condition for flare landing.

3 NSGAII based homing trajectory planning

3.1 Normalization process

In the process of homing trajectory optimization, the significant digit tends to be lost owing to the large orders of magnitude difference of state variables. Normalization can effectively overcome this difficulty and improve the accuracy of calculation. In addition, to keep state variables in the same order of magnitude is often required in the process of homing trajectory optimization. Consequently, new variables are defined as

                             (6)

Then, the motion equations (Eq. (1)) are transformed into

                      (7)

The initial conditions (Eq. (2)) are transformed into

                          (8)

The terminal constraint Eq. (3) are transformed into

                             (9)

The control constraint (Eq. (4)) are transformed into

                      (10)

The three objective function (Eq. (5)) are transformed into

                  (11)

3.2 Parametric method

The essence of homing trajectory planning of parafoil system is a kind of multi-objective optimal control problem. The optimal solution cannot be directly optimized by NSGAII for the search space is a functional space. Consequently, the first thing to do is to transform the optimal control problem into a parameter optimization problem. Shared parametric methods mainly include direct discrete method, multiple parameter interpolation method, function approximation method, and so on. In order to facilitate encoding, decoding and improve the expression ability of control law, non-uniform b-spline is used for fitting the control law in this work. Accordingly, chromosomes of genetic space are formed by control vertices of b-spline basis function. By this method, few dimension control parameters are required for all forms of complex control curves, the encoding, and decoding, as well as the subsequent optimization computation complexity are greatly simplified.

Non-uniform b-spline curve is defined as

                           (12)

where d_i represents control vertices; N_i,k(s) (i={0, 1, …, n}) represents the k-th of b-spline basis function which is decided by s_i (i={0, 1, …, n+k+1}), and can be deduced by De-Boor-Cox recursive formula as

(13)

The support interval of the k-th b-spline basis function is [s_i, s_i+k+1], which contains k+1 node intervals. Namely, the b-spline basis function is related to as much as k+1 nodes, but has nothing to do with the other nodes. The non-uniform vector s=[s₀, s₁, …, s_i+k+1] is decided by Hartley-Judd method as

               (14)

By using real number coding method, chromosomes are formed by d_i, which can be represented as {d₀, d₁, …, d_n}. Generally speaking, the more number the control vertices is, the more precise the control law can be expressed, yet the calculation will be cumbersome, so that the choice of control vertex size needs to balance between the two.

3.3 Improvements of NSGAII

This work proposes an improved NSGAⅡto solve the multi-objective homing trajectory planning problem. The following mainly illustrates the improvements of NSGAII.

3.3.1 Population initialization

The chaos population initialization method is used instead of random population initialization to make the initial population evenly distribute and overcome the premature convergence in calculation process for the particularly inherent randomness and ergodicity of chaos. This paper adopts the logistic chaotic map:

          (15)

where L represents the population size; P_i,N represents random number with 1×N dimension.

3.3.2 Population trimming

For the classical NSGAII, individuals with low rank counter tend to be retained preferentially. Once individuals are trapped in local optimum, it is difficult for them to jump out for all the individuals with high rank counter are pruned away and there is no new individual to update. Thus, in order to maintain the diversity of population and avoid premature phenomenon, a crowding distance based population trimming method is proposed, and the basic idea of the trimming method can be described as follows: a coefficient is introduced such that the number of individuals in the first front that allowed to reserve can be calculated according to its value. While the number of reserved individuals in other fronts allowed to reserve is calculated according to certain formulas. Therefore, the number of reserved individuals in front is the smaller one between the number of individuals allowed to retain and the existing ones. That is to say, for the first front, the number of reserved individuals is directly determined by the coefficient. When the number of reserved individuals is less than the existing ones, the coefficient has a limitative effect on the number of individuals reserved in front. The following work is trimming the number of individuals in front to the number allowed to keep through the tournament selection. As mentioned above, for the same front, the smaller the crowding distance of individual is, the worse its diversity is. Thus, it is more likely to be eliminated.

3.3.3 Constraint handling

According to the practical engineering design requirements, homing trajectory planning of parafoil system usually has specific requests for certain performance indexes, such that, objective functions and constraints have to be redefined according to the requirements. Take one of the practical engineering requirements as an example: the distance error from the landing site to the target is within a circle of specified radius 1 m. In view of this demand, the following can be done: the tolerance △ε is introduced, the objective function is considered a constraint, and then objective function can be transformed into

                    (16)

where c(x) denotes the inequality constraint; △ε can be set according to the practical engineering demand. As for the above example, △ε can be calculated as

                            (17)

Through the above processes, the multi-objective homing trajectory planning problem is transformed into an optimization problem with one constraint and two objective functions. The two objective functions are optimized by NSGAⅡ at the same time, and the constraint can be handled by the static penalty function. The static penalty function is applied to adjust the fitness of individuals as

                 (18)

where is the original fitness; is the fitness after adjustment; P(x) is the penalty function, which is defined as

                              (19)

where k is the coefficient of penalty function and mainly used for excessive or insufficient of punishment and its value is set according to the satisfaction degree of constraints and calculation efficiency, which is usually set to 1.

3.3.4 Optimal solution selecting

In order to avoid possible uncertainties caused by subjective factors of engineers in the process of selecting solutions, a method based on fuzzy set theory is introduced to select the final optimal solution. The main process of this selecting method is that: first, calculate the proportion of the objective function value of a solution to one of the optimization objective. The calculation formula is:

                           (20)

where i indicates the i-th of Pareto optimal solution; j indicates the j-th of objective function; indicates the minimum objective function value of all solution; indicates the maximum objective function value.

And then calculate the comprehensive proportion u_i of each solution of all optimization objectives and its value is limited on [0, 1]. u_i indicates the satisfaction degree of each solution. The larger the value is, the higher the evaluation level of the solution will be. The calculation formula is

                              (21)

where k indicates the number of Pareto optimal solutions in Pareto front; m indicates the number of objective functions.

3.3.5 Algorithm process

The main process of the improved NSGA II can be described as follows: Initially, a chaotic parent population is created. And then the usual binary tournament selection, recombination and mutation operators are used for creating the new offspring population. After that the parent population and its offspring merge into one population. The next step is non-domination sorting, through which, all individual in the population are assigned to different fronts. Then the crowding distances of each individuals in the front are calculated such that the population is trimmed to the same size as the initial population by crowding distance based population trimming method.

4 Numerical simulation

4.1 Parameter settings

The simulation experiment is conducted based on a certain type of parafoil system with the physical parameters shown in Table 3.

Table 3 Design parameters of parafoil system

According to the selected parafoil, basic motion parameters are set as follows under the premise of the bank angle less than 20°. The initial velocity v_s=13.8 m/s, v_z=5 m/s, the initial height h=2000 m and the motor control volume η_max=30. Parameters of the normalized model are set as follows: x₁(0)=0.6, x₂(0)=0.6, x₃(0)=π. The corresponding actual coordinates are x₀=3000 m, y₀=3000 m, ψ₀=π.

In this work, commonly used 5 control vertices and the 2nd of b-spline basis function are used. The improved NSGA II operation parameters are set as follows. The population size is 100; the Pareto optimal solutions number is 30; The number of iterations is 200.

4.2 Simulation results and analysis

In order to comprehensive validate the applicability of the improved NSGA II on homing trajectory planning of parafoil system, three different kinds of schemes are established on the basis of practical engineering demands, as shown in Table 4, where Scheme 3 needs to be specially interpreted. Considering that the parafoil system is dropped near the target, the energy carrying is completely satisfied with all the control needs, such that and are chosen as optimization objective, while the energy consumption is ignored.

Pareto optimal solution sets of homing control laws under three different kinds of schemes obtained by the improved NSGA II are shown in Figs. 1-3.

In the figures above, objective function 1 refers to , objective function 2 refers to and objective function 3 refers to x axis and y axis are respectively corresponding to relevant optimization objects. It can be seen from the figures that all solutions in the Pareto front are well distributed and converged to the Pareto optimal solution.

The final solution of the problem is usually selected according to the following principles: if the y performance requirement of practical problem is high, the solution should be chosen as close as possible to the y axis, so that its x performance is relatively poor, and vice versa. For different design requirements, the improved NSGA II can give all the Pareto optimal solution sets. The final solution can be selected by engineers according to the practical engineering demands which will possibly bring uncertainties due to subjective factors in the selection. Therefore, the method based on fuzzy set theory is applied to select the final optimal solution. Due to the possible uncertainty of constraints handling, the other constraints need to be taken into account when u_i is calculated.

By means of the selecting method, three optimal solutions are selected from each Pareto solution set. By substitution of each optimal solution into the mathematical model of parafoil system, the normalized optimal homing trajectories and corresponding optimal motor control quantity curves are shown in Figs. 4 and 5. It can be seen that homing trajectories and control quantity curves of three schemes with the same initial state and target point are various due to different optimization objects. Homing trajectory curves are relatively smooth, and the gliding period is long. At the gliding stage, the control quantity is relatively small, all of which is intensely well accorded with the homing control characteristics of parafoil system.

Table 4 Planning schemes

Fig. 1 Pareto optimal solution set of Scheme 1

Fig. 2 Pareto optimal solution set of Scheme 2

Fig. 3 Pareto optimal solution set of Scheme 3

Fig. 4 Horizontal homing trajectory of three schemes

Fig. 5 Motor control quantity of three schemes

The total energy consumption, landing distance error and landing direction of the three schemes are shown in Table 5. It can be seen from the data objective functions that have been taken as optimization objectives tend to obtain better optimization precision. Meanwhile, objective functions taken as constraints can meet design requirements with poorer accuracy. The homing results show that the improved NSGA II can extremely well complete homing trajectory planning, and the results fully meet and exceed the practical engineering design requirements. Following the homing trajectory, the parafoil system is able to fly to the target accurately, and implements upwind alignment at the landing time.

Table 5 Homing result

5 Conclusions

1) Variables were normalized into the same order of magnitude so as to improve the accuracy of calculation.

2) Non-uniform b-spline was used for fitting the control law so as to use the algorithm conveniently.

3) Some improvements were made on NSGA II: chaos operator was used for initializing the population to overcome the premature convergence; a crowding distance based method was used in population trimming; the static penalty function with tolerance was introduced to deal with constraints; and the optimal solution was selected in Pareto solution set by method of fuzzy set theory instead of artificial selection so as to obtain the solution with best comprehensive index.

4) Simulation results demonstrate that this method is a kind of effective way for homing trajectory planning with strong robustness and extensive adaptability, and it is more able to reflect the nature of multi-objective optimization. The optimized trajectory makes the parafoil system implement homing in a nearly optimal way. In addition, this method will bring more design choices to practical engineering of homing trajectory planning of parafoil system.

References

[1] LI Chun, LU Zhi-hui, HUANG Wei, SHEN Chao. Guidance navigation & control system for precision fix-point homing parafoil [J]. Journal of Central South University: Science and Technology, 2012, 43(4): 1331-1335. (in Chinese)

[2] ROGERS J, SLEGERS N. Robust parafoil terminal guidance using massively parallel processing [J]. Journal of Guidance, Control, and Dynamics, 2013, 36(5): 1336-1345.

[3] ZHU Er-lin, SUN Qing-lin, TAN Pan-Long, CHEN Zeng-qiang, HE Ying-ping. Modeling of powered parafoil based on Kirchhoff motion equation [J]. Nonlinear Dynamics,2015, 79(1): 617-629.

[4] XIONG Jing. Research on the dynamics and homing project of parafoil system [D]. Changsha: National University of Defense Technology, 2005. (in Chinese)

[5] PEARSON A. Optimal control of a gliding parachute system [R]. Massachusetts: Army Natick Labs, 1972.

[6] GIMADIEVA T. Optimal control of a gliding parachute system [J]. Journal of Mathematical Science, 2001, 103(1): 54-60.

[7] Rademacher B J, Lu P, Strahan A L, Cerimele C J. In-flight trajectory planning and guidance for autonomous parafoils [J]. Journal of Guidance, Control, and Dynamics, 2009, 32(6): 1697-1712.

[8] LIU Zhao, KONG Jian-yi. Path planning of parafoil system based on particle swarm optimization [C]// Computational Intelligence and Natural Computing. Wuhan, China: IEEE, 2009: 450-453.

[9] JIAO Liang, SUN Qing-lin, KANG Xiao-feng. Route planning for parafoil system based on chaotic particle swarm optimization [J]. Complex Systems and Complexity Science, 2012, 9(1): 47-54. (in Chinese)

[10] ZHANG Li-min, GAO Hai-tao, CHEN Zeng-qiang, SUN Qing-lin, ZHANG Xin-hui. Multi-objective global optimal parafoil homing trajectory optimization via Gauss pseudo-spectral method [J]. Nonlinear Dynamics, 2013, 72(1/2): 1-8.

[11] GAO Hai-tao, ZHANG Li-min, SUN Qing-lin, SUN Ming-wei, CHEN zeng-qiang, KANG Xiao-feng. Fault-tolerance design of homing trajectory for parafoil system based on pseudo-spectral method [J]. Control Theory & Applications, 2013, 30(6): 702-708. (in Chinese)

[12] Luders B D, Sugel I, HOW J P. Robust trajectory planning for autonomous parafoils under wind uncertainty [C]// AIAA Conf on Guidance, Navigation and Control and Co-located Conferences. Reston, VA: AIAA, 2013: 1-27.

[13] CLEMINSON J. Path planning for guided parafoils: An alternative dynamic programming formulation [C]// AIAA Aerodynamic Decelerator Systems (ADS) Conference. Reston, VA: AIAA, 2013: 1-21.

[14] BABU A, SUJA V, REDDTY C. Three dimensional trajectory optimization of a homing parafoil [C]// 3rd International Conference on Advances in Control and Optimization of Dynamical Systems. Netherlands: Elsevier, 2014: 847-854.

[15] Tao Jin, SUN Qing-lin, ZHU Er-lin, CHEN Zeng-qiang. Quantum genetic algorithm based homing trajectory planning of parafoil system [C]// 34th Chinese Control Conference. Hangzhou, China: IEEE, 2015: 2523-2528.

[16] DEB K, AGRAWAL S, PRATAP A, MEYARIVAN T. A fast elitist non-dominated sorting genetic algorithm for multi-objective optimization: NSGA-II [J]. Lecture Notes in Computer Science, 2000, 1917: 849-858.

[17] DEB K, PRATAP A. A fast and elitist multi-objective genetic algorithm: NSGA-II [J]. IEEE Transactions on Evolutionary Computation, 2002, 6(2): 182-197.

[18] ABOLFAZL K. Best compromising crashworthiness design of automotive S-rail using TOPSIS and modified NSGA II [J]. Journal of Central South University, 2015, 22(1): 121-133.

[19] LEI Ying-jie, ZHANG Shan-wen. Genetic algorithm toolbox and its application [M]. Xi’an: Xi’an University of Electronic Science and Technology Press, 2014. (in Chinese)

[20] JIAO Liang, SUN Qing-lin, KANG Xiao-feng, CHEN Zeng-qiang, LIU Zhong-xin. Autonomous homing of parafoil and payload system based on ADRC [J]. Control Engineering and Applied Informatics, 2011, 13(3): 25-31.

(Edited by YANG Hua)

Foundation item: Project(61273138) supported by the National Natural Science Foundation of China; Project(14JCZDJC39300) supported by the Key Fund of Tianjin, China

Received date: 2015-08-31; Accepted date: 2015-10-26

Corresponding author: SUN Qing-lin, PhD, Professor; E-mail: sunql@nankai.edu.cn