Controlling underactuated systems is a challenging problem in control engineering. This paper presents a novel constraint-following approach for control design of an underactuated two-wheeled mobile robot (2 WMR), which has two degrees-of-freedom (DOF) to be controlled but only one actuator. The control goal is to drive the 2 WMR to follow a set of constraints, which may be holonomic or nonholonomic constraints. The constraint is considered in a more general form than the previous studies on constraint-following control (hence including a wider range of constraints). No auxiliary variables or pseudo variables are required for the control design. The proposed control only uses physical variables. We show that the proposed control is able to deal with both holonomic and nonholonomic constraints by forcing the constraint-following error to converge to zero, even if the system is not initially on the constraint manifold. Using this control design, we investigate two cases regarding different constraints on the 2 WMR motion, one for a holonomic constraint and the other for a nonholonomic constraint. Simulation results show that the proposed control is able to drive the 2 WMR to follow the constraints in both cases. Furthermore, the standard linear quadratic regulator (LQR) control is applied as a comparison in the simulations, which reflects the advantage of the proposed control.

## Introduction

In recent years, studies on two-wheeled mobile robot (2 WMR) or wheeled inverted pendulum (WIP), especially the motion control for it, have received extensive attention in both academia and industry worldwide [1–4]. This is mainly attributed to its capability to model a class of modern vehicles transporting human or goods safely and efficiently [2], and to demonstrate various control strategies. For instance, SEGWAY, which is a well-known and popular commercial personal transporter, is essentially a 2 WMR; Mu et al. [5], Huang et al. [6], and Pathak et al. [7] applied their proposed control approaches to the WIP for demonstration. From the view point of motion control design, a salient feature of the 2 WMR or WIP is that it is usually underactuated. That is, it has fewer control inputs than the degrees-of-freedom (DOF) to be controlled. This presents a tough problem frequently encountered by researchers, i.e., control design for underactuated mechanical systems [8].

Despite the complexity, the problem of control design for underactuated systems has been among hotspot issues for years. This is due to its wide existence in engineering practice—there are a great number of practical mechanical systems belonging to this community, e.g., underwater vehicles, unicycle robots, humanoid robots, and aircraft. The existing studies on this problem include a variety of methodologies accounting for different issues. In general, these studies can be classified into two types: kinematics approach and dynamics approach. The kinematics approach aims to regulate displacement and/or velocity by controlling other kinematic properties. The control is not force or torque, and examples are Refs. [9,10]. The dynamics approach designs controls are closely related to force or torque, and examples are [11–13]. Some conventional control methodologies also have been applied to underactuated systems [14–16], and most of these studies belong to the dynamics approach. Regarding examples of 2 WMR and WIP, Yang et al. [17] applied neural network-based control approach for an underactuated WIP; Cui et al. [18] investigated adaptive backstepping control for underactuated WIP models; Guo et al. [19] and Xu et al. [20] employed sliding-mode control for an underactuated WMR. Despite the fact that so many studies have been done, the problem of control design for underactuated systems still calls for new developments and more insight. This study targets this problem and aims to design control that falls into the dynamics approach.

In this paper, to design the motion control for an underactuated 2 WMR, we propose a new approach based on constraint-following [21,22]. The 2 WMR consists of two wheels in parallel and an inverse pendulum. The model for the motion of the 2 WMR is formulated as a 2DOF (wheels displacement and pendulum tilting angle) system with one control input (single actuator). The constraints can be holonomic or nonholonomic. The control design is then carried out in two steps. First, a control design without addressing possible initial constraint deviation is proposed. Second, an additional control action is added to deal with possible initial constraint deviation. The control is then obtained as the sum of the controls proposed in the first and second steps. We show that the proposed control is able to force the constraint-following error to converge to zero, even if the system is not initially on the constraint manifold. Following this framework, we investigate two cases with different constraints on the 2 WMR motion. The first one is with a holonomic constraint, which is to render the wheels move forward while keeping the pendulum approximately upright. The second one is with a nonholonomic constraint, which is to render the tilting angle of the pendulum to vary in a specific interval with its midpoint at the upright position. The corresponding control input of the 2 WMR is then obtained based on the constraint by using the proposed control design approach.

The main contributions of this study are threefold. First, we formulate the underactuated 2 WMR motion control problem as constraint-following. The constraints may be holonomic or nonholonomic, thus including a wide range of practical robot needs. Second, we propose a control that can achieve constraint-following for the 2 WMR, even if it is not initially on the constraint manifold. Despite the system being under constraint, no auxiliary variables (such as Lagrange multiplier) or pseudo variables (such as generalized speeds) are needed for the control design. The control only uses physical variables. Third, we investigate two cases with different constraints (one holonomic and one nonholonomic) on the 2 WMR motion and corresponding numerical simulation results are provided to demonstrate the effectiveness of the control.

## Constrained Mechanical System

Here, *t* ∈ **R** is the time, *q* ∈ **R*** ^{n}* is the position vector, $q\u02d9\u2208Rn$ is the velocity vector, $q\xa8\u2208Rn$ is the acceleration vector, and

*τ*∈

**R**

*with*

^{m}*m*≤

*n*is the control input vector. Furthermore, $M(q(t),q\u02d9(t),t),\u2009C(q(t),q\u02d9(t),t),\u2009g(q(t),q\u02d9(t),t)$, and $B(q(t),q\u02d9(t),t)$ stand for the inertia matrix, the matrix of Coriolis/centrifugal terms, the gravitational force, and the input matrix, respectively, and they are of appropriate dimensions. The functions

*M*(⋅),

*C*(⋅),

*g*(⋅), and

*B*(⋅) are continuous. In addition, the position vector

*q*does not need to be the generalized coordinate and can be chosen based on the specifications of the problem.

where $q\u02d9i$ is the *i*th component of $q\u02d9,\u2009and\u2009A\u2032li(\xb7)$ and $c\u2032l(\xb7)$ are both *C*^{1}.

*Remark 1.* The proposed constraints are in the *first-order form*. They may not be integrable and generally may be nonholonomic. Furthermore, the proposed constraints have similar form as those in Ref. [21]. However, they are different. To be specific, the constraints in Ref. [21] do not include the velocity vector $q\u02d9$ as an argument in the functions “$A\u2032li(\xb7)$” and “$c\u2032l(\xb7)$,” but the proposed constraints do. Therefore, the proposed constraints are more general than those in Ref. [21].

*second-order form*, we take derivative of Eq. (2) with respect to

*t*and obtain

and where $c=[c1\u2003c2\u2003\cdots \u2003cm]T$.

Definition 1. *The constraint (9) is consistent if for given A and b, there is at least one solution*$q\xa8$.

Assumption 1. *The constraint (9) is consistent*.

Lemma 1 [23]. *A necessary and sufficient condition for the constraint (9) to be consistent is A ^{+}Ab = b. Here, “+” denotes the Moore–Penrose (MP) generalized inverse*.

The constraint (9) is in fact a very general form. Both holonomic and nonholonomic constraints can be expressed as this form. Task constraints can also be put in this form, e.g., Ref. [24]. Other control problems such as stabilization, trajectory following, and optimality can also be cast into this form by taking derivatives [22].

## Servo Control Design Without Initial Condition Deviation

Definition 2. *For given Φ _{1} and b, the constraint (14) is called consistent if there exists at least one solution*$x\xa8$.

Since there is a one to one correspondence between $x\xa8$ and $q\xa8$, Eq. (9) is consistent if Eq. (14) is consistent.

Definition 3. *The system (13) is called servo constraint controllable with respect to the constraint (14) if there is a control τ _{n} such that the system under this control meets Eq. (14) for all*$(q,q\u02d9,t)\u2208Rn\xd7Rn\xd7R$.

where *τ _{n}* ∈

**R**

*is the unknown. This equation can be viewed as a constraint on*

^{m}*τ*.

_{n}Definition 4. *For given Φ _{1,2} and*$b\xaf$

*, Eq. (15) is consistent if there exists at least one solution τ*.

_{n}Assumption 2. *Equation (15) is consistent*.

*existence*condition for the servo control

*τ*, which can render the nominal system to follow the constraint (9). A direct solution is $\tau n=(\Phi 1\Phi 2)+b\xaf$ such that

_{n}*Subject to Assumption 2, the system (1) is servo constraint controllable with respect to the constraint (14) if and only if*

*for all*$(q,q\u02d9,t)\u2208Rn\xd7Rn\xd7R$

*, and the servo control τ*

_{n}is given by*where S ∈*

*R**$q\u02d9$*

^{m}is an arbitrary vector, possibly dependent on q,*, and t.*

*Proof.*(Sufficiency) Suppose $rank[\Phi 1(q,q\u02d9,t)\Phi 2(q,q\u02d9,t)]\u22651$ for all $(q,q\u02d9,t)\u2208(q,q\u02d9,t)\u2208Rn\xd7Rn\xd7R$, which means the MP inverse of $\Phi 1(q,q\u02d9,t)\Phi 2(q,q\u02d9,t)$ always exists [23] and the control (18) is meaningful. Substituting the control (18) into the system (13) and premultiplying both sides of Eq. (13) by Φ

_{1}yield

Here, a property of MP inverse, *GG*^{+}*G *=* G*, is applied.

_{1}yields

There are two consistency provisions in this section. The first (Assumption 1) ensures the constraint to be meaningful, for otherwise the constraint-following problem is impossible. The second (Assumption 2) ensures the control design to be meaningful, for otherwise there is no way to come up with the control.

*Remark 2.* The servo control *τ _{n}* in Eq. (18) is a

*model-based state feedback control*and is readily applicable. Note that the control is continuous in the state. The condition (17) simply means that the rows of Φ

_{1}and the columns of Φ

_{2}should not be all perpendicular to each other. This allows the control

*τ*to be able to project a specific amount of component in $R(\Phi 1T)$. As shown before, eventually the system needs, after applying the control, a net projection of $\Phi 1+b$ in $R(\Phi 1T)$.

_{n}From Theorem 1, we can find that the control (18) is only to account for the second-order form of the constraint, i.e., Eq. (9) or Eq. (14). That is, the control (18) is only able to guarantee Eq. (9) to be satisfied. If the initial condition of the system, $q(t0)\u2009and\u2009q\u02d9(t0)$ with *t*_{0} the initial time, satisfies the first-order form of the constraint (12), i.e., $A(q(t0),q\u02d9(t0),t0)q\u02d9(t0)=c(q(t0),q\u02d9(t0),t0)$, the control (18) can also guarantee Eq. (12) to be satisfied, since Eq. (12) can be obtained by computing the integral for Eq. (9) with respect to *t*. However, it is likely that the initial condition does not satisfy the constraint (12) in practice. Notice that the constraint (12) is the purpose of control design, it is necessary to add an additional control for dealing with possible initial condition deviation so that the constraint (12) can be satisfied, which will be discussed in the following section.

## Dealing With Initial Condition Deviation

Since the constraint (12) is the control purpose, *β* in Eq. (21) can be viewed as a directly measure of the system performance. This will be used in later development.

*Remark 3.* Note that the matrix *ADB*(*ADB*)^{T} is always positive semidefinite for all $q\u2208Rn,\u2009q\u02d9\u2208Rn,\u2009t\u2208R$, hence all eigenvalues are non-negative. What this assumption adds is that the minimum eigenvalue will not equal zero and is always a finite distance away from 0 for all $q\u2208Rn,\u2009q\u02d9\u2208Rn,\u2009t\u2208R$ (therefore not infinitesimally close to 0).

where *κ* > 0 is a constant.

Theorem 2. *Subject to Assumptions 2 and 3, consider the system (1). The control (25) renders the following performance:*

- (i)
Uniform stability: For each

*ζ*> 0, there exists*ξ*> 0, such that if*β*(⋅) is any solution with $\Vert \beta (t0)\Vert <\xi $, then $\Vert \beta (t)\Vert <\zeta $ for all*t*>*t*_{0}. - (ii)

*V*with respect to

*t*yields

Since the upper bound of the derivative of *V* is nonpositive and equals 0 only if *β* = 0, we have the uniform stability as well as the convergence of *β* to 0 as *t* → *∞*. ◻

## Dynamical Model of an Underactuated Two-Wheeled Mobile Robot

Up to now, many dynamical models for two-wheeled mobile robots have been proposed. This paper adopts the one in Ref. [20], which takes frictions and input coupling into consideration. Figure 1 depicts the schematic of the 2 WMR. The wheels move along the road surface, and their displacement and velocity are, respectively, denoted by *x* and $x\u02d9$ with rightward as the positive direction. The tilting angle and the angular velocity of the pendulum are, respectively, denoted by *θ* and $\theta \u02d9$, with the upright position as origin and clockwise rotation as the positive direction. *φ* is the slope angle of the road, and for the flat road, *φ* = 0. *f _{r}* is the friction between the wheels and the road.

*f*is the joint friction acting on both the wheels and the pendulum as

_{j}*f*and −

_{j}*f*, respectively.

_{j}*u*is the control torque generated by the motor which directly acts on the wheels with clockwise rotation as the positive direction. As the motor is directly mounted on the pendulum, there is a reaction torque −

*u*directly applied to the pendulum.

*m*,

_{w}*I*, and

_{w}*r*are the mass, rotation inertia, and radius of the wheels, respectively,

*m*and

_{p}*I*are the mass and rotation inertia of the pendulum, respectively,

_{p}*l*is the distance between center of gravity of the pendulum and the center of the wheel, and

*g*

_{0}is the acceleration of gravity. The frictions are modeled as a combination of viscous friction and Coulomb friction. That is, $fr=fvx\u02d9+fcsgnx\u02d9\u2002and\u2002fj=fjv\theta \u02d9+fjcsgn\theta \u02d9$, where

*f*and

_{v}*f*are viscous-friction constants,

_{jv}*f*and

_{c}*f*are Coulomb-friction constants, and sgn(⋅) represents the signum function. Substituting the expressions of the frictions into Eqs. (32) and (33) yields

_{jc}## Control Design for the Two-Wheeled Mobile Robot

### Constraint on the Two-Wheeled Mobile Robot Motion.

The proposed control design in Secs. 3 and 4 requires the number of the constraint to be equal to that of the control input. Since there is only one available control input, it is only able to impose one constraint on the 2 WMR motion. Generally, there are many choices of the constraint that can be expressed as Eqs. (12) and (9), either holonomic or nonholonomic. For well demonstration of the proposed control, we will consider two cases regarding different constraints—a holonomic one and a nonholonomic one.

#### Case 1: Holonomic Constraint.

where *θ _{d}* > 0 is a constant. This constraint can render the wheels to move forward while keeping the pendulum approximately upright if we impose some restrictions on

*θ*. The arguments are as follows.

_{d}*θ*= 0,

*θ*should be small (hence be close to “0”) so that the pendulum can be “approximately upright.” Second, taking derivative of Eq. (37) with respect to

_{d}*t*twice yields

*θ*is small (close to zero) and the absolute value of the slope angle of road, $|\phi |$, is usually much smaller than

_{d}*π*/2, it is reasonable to assume that $\u2212\pi /2\u2264\theta d+\phi \u2264\pi /2$. Therefore, we have $a\u20322\u22650$. Now we impose the following restriction on

*θ*:

_{d}*r*,

*m*,

_{p}*m*, and

_{w}*l*are the designer's discretion. Consequently, we have

*r*,

*a*

_{1}, $fv>0,\u2009a\u20322\u22650$, and $mplg0\u2009sin\u2009\theta d\u2212rfc\u2212r\u2009sin\u2009\phi (mp+mw)g0>0$, we can obtain from Eq. (43) that $x\u02d9$ will converge to a positive value, which is

Since $mplg0\u2009sin\u2009\theta d+rfc\u2212r\u2009sin\u2009\phi (mp+mw)g0>0$, we can obtain from Eq. (45) that $x\u02d9$ will converge to a positive value. When $x\u02d9$ becomes positive, Eq. (40) becomes Eq. (43), and consequently, $x\u02d9$ will converge to *v _{d}*. In a conclusion, under the constraint (37),

*θ*will converge to

*θ*, $x\u02d9$ will converge to a positive value, which is

_{d}*v*. That is, the 2 WMR will move forward while keeping the pendulum approximately upright under the constraint (37). Furthermore, since

_{d}*v*depends on

_{d}*θ*, we can adjust the moving speed

_{d}*v*by choosing different

_{d}*θ*.

_{d}Note that the matrix *A* is of full rank. Therefore, $A+Ab=(ATA)\u22121AT\ufe38A+Ab=b$, which verifies Assumption 1 according to Lemma 1. That is, the constraint (46) is consistent.

#### Case 2: Nonholonomic Constraint.

where *k*_{0} > 0 and *k*_{1} ≥ 1 are constants. This constraint implies that the moving direction of the wheel will be consistent with the tilting direction of the pendulum. Furthermore, as the domain of $tan(\xb7)$ is discontinuous at $k\u0303\pi +\pi 2\u2009(k\u0303\u2208Z)$, we can choose $k1\theta (t0)\u2208(\u2212\pi 2,\pi 2)$ so that $k1\theta (t)\u2208(\u2212\pi 2,\pi 2)$ for all *t* > *t*_{0}. Therefore, this constraint can render the tilting angle of the pendulum *θ* to vary in the interval $(\u2212\pi 2k1,\pi 2k1)$, whose midpoint is at the upright position (*θ* = 0). This will be reflected in the simulation results.

Note that the matrix *A* is of full rank, which can verify Assumption 1 by invoking the same argument in the case 1. Therefore, the constraint (48) is consistent.

### The Control Torque.

Since the dynamical model of the 2 WMR can be put in the form of Eq. (1) and the constraints can be put in the form of Eqs. (12) and (9), the robot motion control can be designed according to Secs. 3–5. If there is no initial condition deviation, the control torque, *u*, can be derived as *τ _{n}* in Eq. (18); if initial condition deviation exists,

*u*can be derived as

*τ*in Eq. (25).

### Design Procedure.

Figure 2 shows a flowchart to sum the control design procedure for the underactuated 2 WMR.

## Numerical Simulations

This section presents the numerical simulation results on the 2 WMR to verify the proposed control design approach. As a comparison, the standard linear quadratic regulator (LQR) (linear quadratic regulator) control is also applied to the 2 WMR, and corresponding simulation results are presented. The most common robustness measures attributed to the LQR are a one-half gain reduction in any input channel, an infinite gain amplification in any input channel, or a phase error of ±60 deg in any input channel. The calculation for the standard LQR control gain is based on linear system model. Therefore, to obtain the LQR control gain, we first apply standard (Jacobian) linearization on the dynamical model of the 2 WMR, i.e., Eqs. (32) and (33), and the resulting linear model is then used for the LQR design by invoking the standard procedure of calculating the LQR control gain (selecting the weighting matrices, solving the Riccati equation, etc.). All the simulations are implemented by using the ode15i algorithm in matlab, and the tolerance is chosen to be 10^{−8}.

where $a2=mpl\u2009cos\u2009\theta $. The parameters of the 2 WMR are chosen as follows: $mw=4\u2009kg,\u2009Iw=0.08\u2009kg\xb7m2,\u2009mp=8\u2009kg,\u2009Ip=0.2\u2009kg\xb7m2,\u2009r=0.2\u2009m,\u2009l=0.25\u2009m,\u2009fv=0.5$, *f _{c}* = 1,

*f*= 0.2,

_{jv}*f*= 0.3, and

_{jc}*g*

_{0}= 9.81 m/s

^{2}. Here, the considered 2 WMR is not meant to emulate a SEGWAY. That is, it is unnecessary to consider the 2 WMR with the same size or mass as a SEGWAY. Actually, the 2 WMR model is suitable for various types of vehicles, from outdoor to indoor transportation. The parameters of the 2 WMR can take different values based on different engineering applications. The initial conditions are $x\u02d9(0)=0.1\u2009m/s,\u2009x(0)=0\u2009m,\u2009\theta \u02d9(0)=0\u2009rad/s,\u2009and\u2009\theta (0)=0.08\u2009rad$. Note that the initial conditions do not satisfy the constraints (46) and (48) in the two cases. As discussed in Sec. 6.2, we should apply the control torque

*u*=

*τ*+

_{n}*τ*to the 2 WMR.

_{a}### Simulations for Case 1.

*r*and

*a*

_{1}are positive constants, hence we can obtain

_{1}Φ

_{2}is always a 1 × 1 invertible matrix, indicating that Eq. (15) is consistent and Assumption 2 is verified. Furthermore, Eqs. (51) and (52) imply

*κ*= 0.1. The desired tilting angle

*θ*is small so that the pendulum can be approximately upright. Furthermore, we have $sin\u2009\theta d=0.0998>(rfc+r\u2009sin\u2009\phi (mp+mw)g0/mplg0)=0.0102$, which satisfies the restriction (41). Figure 3 shows the time histories of $\Vert \beta \Vert $ under no control, the LQR control, and proposed control. The simulations under no control and the LQR control suffer interruptions due to the miracle of overflow. We can find that $\Vert \beta \Vert $ histories under no control and the LQR control are far away from zero, while the $\Vert \beta \Vert $ history under the proposed control converges to zero after some time, which echoes Theorem 2. Figure 4 shows the

_{d}*θ*histories. Apparently, only the

*θ*history under the proposed control converges to the desired titling angle

*θ*, the other two histories do not. Figure 5 shows the $x\u02d9$ histories. It can be found that the speed $x\u02d9$ under the proposed control increases in the beginning and finally converges to

_{d}*v*, which is calculated from Eq. (44) as

_{d}*v*= 19.5873 m/s > 0. This verifies the discussions for the constraint in Sec. 6.1.1. Figure 6 shows the histories of control torque

_{d}*u*. The LQR control input appears to be very large (at the level of 10

^{4}N⋅m) at around

*t*=

*6 s, while the proposed control input keeps small all the time. The proposed control torque can be generated by using an appropriate motor. It can be found that the maximum value of the control torque is less than 3 N⋅m. Therefore, the proposed control torque can be generated by a motor whose rated torque is around or larger than 3 N⋅m.*

Using the proposed control, we further investigate the relation between the speed $x\u02d9$ and the desired titling angle *θ _{d}*. Figure 7 shows three different $x\u02d9$ histories under three different desired tilting angles, i.e.,

*θ*= 0.1 rad,

_{d}*θ*= 0.15 rad, and

_{d}*θ*= 0.2 rad. It can be found that the larger

_{d}*θ*is, the faster $x\u02d9$ increases and the larger $x\u02d9$ is. This is actually reflected by Eqs. (43) and (44), which implies that the acceleration and speed of the 2 WMR can be tuned by choosing different desired tilting angles.

_{d}### Simulations for Case 2.

We choose three groups of constraint parameters for demonstration: *k*_{0} = 1, *k*_{1} = 5; *k*_{0} = 1, *k*_{1} = 10; and *k*_{0} = 1, *k*_{1} = 15. It can be noted that the three groups have the same *k*_{0} and different *k*_{1}. The control parameter is chosen as *κ* = 1. Figure 8 shows the time histories of $\Vert \beta \Vert $ for the three groups of constraint parameters. The three $\Vert \beta \Vert $ histories all converge to zero after some time, which echoes Theorem 2. Figures 9 and 10, respectively, show the *θ* and $x\u02d9$ histories of the three groups of constraint parameters. We can find that (1) all *θ* and $x\u02d9$ histories are strictly positive, indicating that the moving direction of the wheel is consistent with the tilting direction of the pendulum and (2) all *θ* histories vary in the range of $(\u2212\pi 2k1,\pi 2k1)$. This result verifies the discussions for the constraint in Sec. 6.1.2. Figure 11 shows the time histories of control torque *u*, which are all mild and can be generated by a motor whose rated torque is around or larger than 6 N⋅m.

## Conclusions

Controlling underactuated mechanical systems is often a very challenging problem. Past researches show many endeavors on applying conventional control methodologies (e.g., sliding model control). This paper, on the other hand, proposes a constraint-following approach to design the control. The constraints may be holonomic or nonholonomic, thus including a wide range of practical system needs. The control is designed in two steps. First, a control without considering initial condition deviation is designed. This sets the base for the second step, which deals with possible initial condition deviation. The control finally obtained is a model-based state feedback control, which is readily applicable. One salient feature is that no auxiliary variables (such as Lagrange multiplier) or pseudo variables (such as generalized speeds) are needed for the control. The control only uses physical variables. We believe that the proposed control design is the simplest and the approach may open a new avenue for controlling underactuated mechanical systems.

Using the proposed control design approach, we investigate two cases regarding different constraints on the 2 WMR motion: one for holonomic constraint and the other for nonholonomic constraint. Both constraints have specific physical meanings. Simulation results on the 2 WMR provided verify the effectiveness of the proposed approach.

## Funding Data

National Natural Science Foundation of China (Grant No. 11572121; Funder ID: 10.13039/501100001809)

Independent Research Projects of State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body in Hunan University (Grant No. 71375004; Funder ID: 10.13039/501100003824)

China Scholarship Council (Grant No. 201606130100; Funder ID: 10.13039/501100004543)

Hunan Provincial Innovation Foundation for Postgraduate (Grant No. CX2016B081; Funder ID: 10.13039/501100010083).