Dynamics and Switching Control of a Class of Underactuated Mechanical Systems with Variant Constraints

: Locomotion systems with variant constraints are familiar in real world applications, but the dynamics and control issues of variant constraint systems have not been sufﬁciently discussed to date. From the viewpoint of Lagrange–d’Alembert equations with additional variable constraints, this paper investigates the modeling approaches of a class of hybrid dynamical systems (HDS) with instantaneously variant constraints and the switching control techniques of stabilizing the HDS to given periodic orbits. It is shown that under certain conditions there possibly exist zero impact periodic orbits in the HDS, and the HDS can be stabilized to the period-one orbits by a linear controller with only partial state feedback, even though the HDS are generally underactuated nonholonomic systems. As an example, a one-legged planar hopping robot is employed to demonstrate the main results of modeling and control of a class of HDS.


Introduction
Variant constraint systems, such as structurally reconfigurable systems, soft body locomotion systems, and contact or collision systems and so on, are familiar in real world applications. The variant constraint systems (VCS) are generally heterogeneous multi-mode systems and contain state jumps in the dynamics of the systems. Thus, the VCS is essentially a class of hybrid dynamical systems (HDS) [1,2].
Although the investigations about the HDS have been ongoing for more than thirty years [3], the dynamics modeling of general HDS has not been thoroughly discussed so far. For the VCS, a reasonable model can be used to analyze the dynamical behaviors of the systems, which can provide useful information for improving the mechanisms/structure designs and simultaneously reducing the complexity of synthesizing the control laws. Unfortunately, to date, a systemic approach of modeling the general VCS has not been presented. As shown in [1], many examples of HDS were presented, but the issues of modeling hybrid systems were discussed case by case. Hyon and Emura [4] briefly presented the modeling approach of the hybrid dynamics of a passive one-legged running robot, and did not mention the theoretical basis of the approach. In a monograph [5], Sadati et al. introduced the impact model of several kinds of dynamical legged robots in brief and the main contents of the monograph are motion planning and hybrid control issues of the HDS. In literature [6,7], Hurmuzlu et al. addressed the modeling and control approaches of the HDS with the background of the biped robots. Following the modeling approach presented by Hurmuzlu [6], Grizzle et al. [8,9] stablized the HDS of biped robots with the finite time controllers. From a viewpoint of nonlinear where L(q,q) = T(q,q) − V(q) is the Lagrangian, T(q,q) is the kinetic energy, V(q) is the potential energy of the system, q ∈ n is the generalized coordinates, λ ∈ m is a vector of Lagrange multipliers, W(q) ∈ n×m is a functional matrix, B ∈ n×s is the input matrix, and u ∈ s is the control input of the system. The term W(q)λ in the right side of Equation (1) denotes the generalized forces caused by external constraints, which are usually introduced by certain interaction between the mechanical systems and the external environment, such as sliding, rolling, elastic or inelastic collision and so on, thus the external constraints are variable. Under the case W(q)λ = 0, the corresponding system (1) is called under the least constraint mode. Otherwise, for the case W(q)λ = 0, the system is called under the full constraint mode. When W(q)λ = 0, as that systematically discussed in [11,24], there necessarily exist m constraint equations If the constraint Equation (2) is integrable, then there exist m algebra equations h(q) = (h 1 , h 2 , . . . , h m ) T = 0 that satisfy W(q) T = ∂h(q) ∂q . (3) Then, system (1) is called a holonomic system. Otherwise, system (1) is a nonholonomic system if the differential constraints Equation (2) are non-integrable.
In order to analyze the state jumps before and after the constraint variation, define t + vc and t − vc to be the time just after and just before the constraint variations, respectively, ∆t vc = t + vc − t − vc is the interval of constraint variation, and defineλ = t + vc t − vc λdt (4) to be the vector of the impulse caused by the constraint variation. In this paper, it is also supposed that the constraint variation is instantaneously completed, that is, ∆t vc ≈ 0, and dq(t) dt ∈ F(t, q(t)), |F(t, q(t))| ≤ δ for all t ∈ t − vc t + vc , where δ ≥ 0 is a positive constant, dq(t) dt ∈ F(t, q(t)) is called the differential inclusions [21]. According to the existence theorem of solutions of discontinuous systems (referring to Theorem 1, chapter 2 of [21]), the solution q(t) of the differential inclusions Equation (5) is absolutely continuous, that is, at the instant t vc , the coordinate variables q(t vc ) do not change For system (1) with the variant constraints Equation (2), the impulse Equation (4) can be analytically calculated due to the following result.

Theorem 1. (Impulse calculation for general systems with variant constraints)
For system (1) with variant external constraints Equation (2), if ∆t vc ≈ 0, then the impulse Equation (4) caused by the constraint variation is explicitly given as for adding constraints, andλ for reducing constraints.
Proof. It is intuition that the following equation is satisfied: In addition, considering assumption Equation (5) and using condition Equation (6), we can conclude that ∂L(q,q) ∂q = 0 and ∂W(q) ∂q = 0 for all t ∈ ∆t vc since the Lagrangian L(q) and the constraints W(q) are invariant during the instant t ∈ ∆t vc . By integrating system (1) on every interval ∆t vc [6,7], we have where Combining with Equations (11) and (12), it follows thaṫ From Equation (2), we have for adding constraints, and for reducing constraints. Substituting Equation (13) into Equation (14), it follows that Then, Equation (7) can be obtained. Substituting Equation (13) into Equation (15), it follows that Then, Equation (8) can be obtained. This completes the proof.
For the underactuated systems considered in this paper, we partition the generalized coordinates into q = (q x , q s ) ∈ n−l × l , of which q x is called the external variables and the external variables are assumed to be passive, and q s is called the shape variables and the shape variables are assumed to be actuated. The shape variables are defined to be the variables that are shown in the inertia matrix M(q s ) of the system. Otherwise, the coordinates q x that are not shown in M(q s ) are defined to be the external variables. Based on the definitions suggested by Reza Olfati-Saber [25], system (1) can be partitioned as where u ∈ s , W x (q) ∈ (n−s)×m , W s (q) ∈ s×m , and both of the equations ∂T ∂q x = 0 and ∂V ∂q = 0 are utilized in Equation (18). In addition, an assumption B = (0, I l ) ∈ n×s is also considered in Equation (18) without losing any generality. For the underactuated system (18) with variant constraints Equation (2), we propose the following corollary based on Theorem 1.

Corollary 1. (Impulse calculation for underactuated systems with variant constraints)
If all the external variables of the underactuated mechanical system (18) are passive, all of the shape variables are actuated, and the constraint variations are instantaneously completed, i.e., ∆t vc ≈ 0, then, for the underactuated system (18), the impulse defined by Equation (4) can be calculated as for adding constraints, andλ for reducing constraints. (4) is an important parameter for analyzing the dynamical behaviors or designing a closed-loop controller for a variant constraint system, which is generally a hybrid dynamical system. Since the impulseλ could not be directly calculated by its definition Equation (4) due to the uncertainties of the time interval ∆t vc = t + vc − t − vc in practice, Equations (19) and (20) provide a feasible way to explicitly calculate the impulseλ caused by the constraint variations for underactuated systems, or Equations (7) and (8) for general dynamical systems.
For many applications, energy efficiency is an important measurement index, and it is necessary to analyze energy changes caused by the constraint variations. The changes of the kinetic energy of the underactuated system (21), (22), (23) or (24) due to the state jumps can be written as By substituting Equation (13) into Equation (25) and applying Equations (14) and (15), it is easy to show that From Equation (26), a useful result for the problems of motion planning and controller design of the hybrid dynamical system (21), (22), (23) or (24) can be presented as follows.
Proof. Supposeλ(t vc ) = 0 for all t ∈ t − vc t + vc , referring to Equation (26), we can conclude that ∆T(q,q) = 0 since the inertia matrix M(q s ) is a positive definite matrix and W(q) is not equal to 0 for a variant constraint system. Accordingly, if the kinetic energy ∆T(q,q) = 0 on ∀t ∈ t − vc t + vc , then there must be ∆T(q,q) < 0 since the trajectory x(t) is governed by the passive dynamics of systems (21), (22), (23) or (24). At the moment, it is necessary that lim Consequently, because of the inertia matrix M(q s ) = 0 for all t ∈ [0, ∞), we can conclude that lim This contradicts the given condition lim Thus, the proof is completed.

Remark 2.
It is worth pointing out that the hybrid dynamical system (21), (22), (23) or (24) is obtained based on the assumption ∆t vc ≈ 0. Otherwise, the dynamical systems should be expressed in different forms. However, modeling the dynamics of the slow varying constraint systems is still an open problem so far.

Remark 3.
On a passive dynamics trajectory x(t), there is no energy consumption because of u(t) = 0. Thus, a nontrivial periodical trajectory x(t) governed by the passive dynamics of a hybrid dynamical system (21), (22), (23) or (24) is useful for all mechanical systems if the time-varying trajectory x(t) is physically meaningful.

Control of Underactuated Systems with Variant Constraints
In this section, we discuss the control issues of the underactuated hybrid dynamical system (21), (22), (23) or (24). Due to the Brockett's theorem [12], which points out that a nonholonomic system can not be stabilized by any smooth time invariant pure state feedback, now it is well known that controller design for an underactuated system with a single mode is not even a simple task. Thereafter, to date, almost all controllers presented in literature for nonholonomic systems adopt time-varying feedback or non-smooth feedback. The underactuated systems (21), (22), (23) or (24) are generally some kinds of nonholonomic systems. This is due to the fact that there are some passive degrees of freedom (DOF) in the systems. Regardless of the control input, we have d dt from Equation (21), and d dt from Equation (22), and the differential Equations (28) and (29) are not integrable generally (except for some special circumstances such as the external variants q x happen to be cyclic coordinates). For the underactuated system (21), (22), (23) or (24), from the viewpoint of control, Equations (28) and (29) are additional differential constraints with respect to the inputs if all state variables of the underactuated systems (21), (22), (23) or (24) are simultaneously stabilized. However, from the viewpoint of mechanics, an underactuated system with partial passive DOF usually satisfies the Lagrange equations without the need for using any external differential constraints [24]. This is different from a traditional nonholonomic system such as a system with rolling constraints that requires the use of Lagrangian multipliers in establishing the dynamics [11]. We notice that a traditional nonholonomic system commonly involves external differential constraints. However, an underactuated system brings about nonholonomic constraints from its own dynamics, such as the underactuated systems that satisfy the law of angular momentum conservation. Therefore, it is more reasonable that a traditional nonholonomic system is renamed as the exogenous nonholonomic systems while an underactuated nonholonomic system with partial passive DOF is named as the endogenous nonholonomic systems. Before discussing the control issues of the systems (21), (22), (23) or (24), it is necessary to point out that the Lagrange multipliers λ in Equation (22) can commonly be eliminated by properly using the dynamics under the least constraint mode and the constraint Equation (2). Then, the state space dimensions of the system (22) under full constraint mode will reduce. Nevertheless, after eliminating the Lagrange multipliers λ, the dynamics of full constraint mode can also be expressed as the form of Equation (21), as it should be, and the dimensions of the state space of the continuous time subsystems in double modes are usually different. Now, taking the control problems of the hybrid systems into consideration, some basic assumptions are given as follows: Basic assumptions: (A1) The actuated DOF s of the underactuated system (21), (22), (23) or (24) in continuous-time modes is not less than the passive DOF n − s; (A2) All the passive subsystems (q x ,q x ) and the actuated subsystems (q s ,q s ) of the hybrid system i.e., ∆t vc ≈ 0.
Owing to the assumptions (A1) and (A2), without loss of generality, both the passive subsystems and actuated subsystems of Equations (21) and (22) can be changed to a linear system with proper dimension by certain globally diffeomorphic coordinate transformations and input changes [26], and the corresponding equivalent system can be expressed as the following switched linear systems [22] x where σ(t) ∈ H, H = {1, 2, . . . , N m } is a piecewise constant function of time, called a switching signal, and N m is the total number of the continuous time subsystems of the hybrid dynamical system (21), (22), (23) or (24). For the passive subsystems of the hybrid dynamical system (21), (22), (23) or (24), the state variables in Equation (30) are given by x σ = col(q x ,q x ) and the equivalent inputs are v σ =q x , while, for the actuated subsystems of the hybrid dynamical system (21), (22), (23) or (24), the state variables and the inputs are given by x σ = col(q s ,q s ) and v σ =q s , respectively. For ∀σ(t) ∈ H, the two matrices A σ and B σ are constant and the pair (A σ , B σ ) is controllable. We consider the trajectory tracking problems for the hybrid system (21), (22), (23) or (24).
to be the instant of the state jumps, where τ * is a constant, and a = N σ ∑ i=1 ∆t i ; then, based on the basic assumptions, a switching controller for the hybrid systems (21), (22), (23) or (24) can be presented as follows.
Theorem 3. If the hybrid systems (21), (22), (23) or (24) satisfy the basic assumptions, and, for a given state trajectory there exists a set of piecewise controllers v σ(∆t i ) (t), σ(∆t i∈S ) ∈ H so that (A4) on the time intervals t ∈ ∆t i∈Ω⊂S and Ω = ∅, the errors of the closed-loop subsystem e s = x s (t) − x d s (t) controlled by the inputs v σ(∆t i ) (t), t ∈ ∆t i∈Ω⊂S are globally exponentially stable and meanwhile the errors of the open-loop subsystem e  (21), (22), (23) or (24) are slow switching systems, that is, the constant τ * is sufficiently large.
Proof. For the linear systems (30), the corresponding error systems can be rewritten aṡ where Then, by the feedbackṽ σ = −K σ e σ , the error systems (31) can be written asė where the matricesĀ σ = A σ − K σ are Hurwitz matrices. Thus, for all t ∈ ∆t i∈S , the error systems of the closed-loop subsystems satisfy where ρ i > 0 are positive constants that are determined by the characteristic values of the matricesĀ σ ; meanwhile, according to (A4) and (A5), the errors of the open loop subsystems satisfy where δ i > 0 are constants that are caused by the state jumps of the hybrid systems (21), (22), (23) or (24), and the constants δ i are finite since the matrix M diag (q s ) is always positive definite and the norm W(q) has an upper bound for all smooth constraints Equation (2). On the other hand, due to Ω = ∅ and S − Ω = ∅, on a subsequent interval ∆t j∈S,j =i , the state errors δ i caused by the i-th state jump can be eliminated because of the assumption (A6), which supposes that τ * is sufficiently large. Then, for the given target trajectory x d (t) on t ∈ [0, a], we can conclude that there exist constants 0 <ρ i < ρ i , i ∈ S so that This completes the proof.

Remark 4.
Theorem 3 reveals that a class of the underactuated systems with variant constraints can track certain time varying trajectory by switched linear controllers with only partial state feedback in a time-sharing manner, even though the underactuated systems under consideration are generally nonholonomic systems that are possibly caused by either external non-integrable differential constraints or internal dynamics of the passive subsystems (28) and (29). Thus, the presented switched linear controllers are essentially a class of discontinuous nonlinear feedback controllers [27]. It is interesting that the HDS (21), (22), (23) or (24) could be stabilized to a time-varying trajectory by the presented switched linear controller, but the presented controllers could not stabilize the systems (21), (22), (23) or (24) to time invariant equilibrium points since the full state variables should be stabilized simultaneously.
For many applications in practice, energy-efficiency is an important index. On the other hand, from a viewpoint of the closed loop systems, the larger energy consumption in steady motions of a controlled system often means the less of stability margin of the system. Thus, it is useful to search for an energy efficient trajectory of the hybrid systems (21), (22), (23) or (24) so that the stabilizing conditions of the switched linear controllers presented in Theorem 3 could be further relaxed. To this end, we present the following statement. Then, systems (21), (22), (23) or (24) can be globally exponentially stabilized to x d (t) by the switched controllers v σ(∆t i ) (t) according to the switching sequence (σ(∆t 1 ), σ(∆t 2 ), . . . , σ(∆t N σ )).
Proof. Referring to systems (23) and (24) and considering the condition (B1), we can see that the state jumps are zero. This is possible since the solutions (q,q) of the following equations and generally exist and not unique. Without loss of generality, suppose the closed-loop systems on the intervals ∆t i∈Ω⊂S happen to be the passive subsystems. According to the condition (B2), and referring to Equation (33), the errors of the closed-loop subsystems satisfy where ρ i > 0 are positive constants that can be determined by the characteristic values of the matrices A σ ; meanwhile, the errors of the open loop subsystems have relationships where γ i > 0 are also positive constants. On the other intervals, ∆t i∈S−Ω , the closed-loop subsystems should be the actuated subsystems, and the errors of the closed-loop subsystems satisfy where α i > 0 are constants; meanwhile, the errors of the open loop subsystems have relationships where β i > 0 are also constants. Note that the four groups' constants ρ i , γ i , α i and β i satisfy the following relationships because of the given condition ∑ i∈Ω and ∑ i∈S−Ω Therefore, on the interval [0, a], the errors of systems (21), (22), (23) or (24) can be estimated by where and From Equation (43), it is directly shown that ρ 0 < 0, and γ 0 < 0. This completes the proof.
Remark 5. Theorem 4 shows that the stabilizing conditions presented in Theorem 3 can be relaxed if the motions of the controlled system satisfy the zero impulse condition Equation (27) [28]. In Theorem 4, the open loop subsystems are permitted to be unstable while they satisfy certain conditions. More importantly, by combining Theorems 2 and 4, it is not difficult to find that the hybrid dynamical system (21), (22), (23) or (24) could be stabilized to a smooth periodic orbit that is governed by the nature dynamics of the system while the periodic orbit also satisfies the zero impulse condition Equation (27).

The Hybrid Dynamics of a Planar One-Legged Hopping Robot
In order to demonstrate the applications of the main results presented in the preceding sections, a planar hopping robot system [4], as shown in Figure 1, is employed to achieve this goal. The hopping robot has a single actuated linear spring leg and a rotational actuated body, and a torsion spring is installed in the hip joint. Suppose the mass of the leg and the body are m l and m b , respectively, while the foot is massless. The moment of inertia of the leg and the body are J l and J b , respectively. The length of the leg is l, and its nature length is l 0 . The spring stiffness of the leg is k l and the stiffness of the hip spring is k h . For the purpose of simplicity, the center of mass (CM) of the body and the leg are designed just on the hip joint, and the mass of the springs is supposed to be negligible. The generalized coordinates of the hopping system are defined to be (x, z, θ, ϕ, l), where (x, z) denotes the position of the foot tip, θ represents the angle between the leg and the horizontal plane, ϕ represents the angle between the body and the leg, and the positive direction of all angular variable is defined to be counter clockwise.
The dynamics of the hopping robot in flight phase can be expressed as where τ f is the input torque of the hip joint, and the position of the CM of the robot is given as Note that, in flight phase, the length of the leg is l = l 0 since the leg can not be controlled in flight phase. Nevertheless, in stance phase, Equations (48) are the holonomic constraints of the hopping robot. Thus, the Jacobian matrix W(q) of Equation (2) can be written as follows by Equation (3): and the constraint forces are Thus, the dynamics of the hopping robot in stance phase can be written as By the constraint Equation (48) and the first two equations of Equation (51), the constraint forces (λ x , λ z ) can be resolved, and then the dynamics Equation (51) can be simplified to We can verify that the dynamics Equation (52) are the same as that derived directly from the Euler-Lagrange equations. In order to get the equations of the state jumps for the hopping robot, from Equation (11), we have where λ x = −p sin (θ) , since the impulse p is perpendicular to the longitudinal axis of the leg, while, along the axis of the leg, the impulse can be absorbed by the springs of the robot. Substituting Equation (54) into Equation (53), it follows that Then, during the constraint variation interval ∆t vc = t + vc − t − vc , the changes of the speed of the generalized coordinates can be written from Equation (55) as where θ(t ± vc ) denotes θ(t + vc ) = θ(t − vc ). On the other hand, from the constraint Equation (48), it is easy to show From Equation (57), it can be deduced thaṫ Substituting Equation (56) into Equation (58) and doing some simple calculations, the impulse can be calculated as where Then, substituting Equation (59) into Equation (56), the equations of state jumps caused by the constraint variations can be given as andl(t + vc ) =ẋ cm (t + vc ) cos θ(t ± vc ) +ż cm (t + vc ) sin θ(t ± vc ) for adding constraints due to Equation (57), andl(t + vc ) = 0 for reducing constraints. Thus, the hybrid dynamics of the hopping robot are composed of Equations (47), (52) and (61).

Motion Planning for the One-Legged Hopping Robots in Steady Locomotion
By introducing springs into the hopping robots, as shown in [29], it is possible for the robots to realize steady hop motions without using large control inputs, while the steady motions are primarily governed by the nature dynamics of the dynamical systems. As seen in the last subsection, the hopping robots are hybrid dynamical systems and the variant constraints of a locomotion system generally result in large state jumps that can easily destroy the stability of the controlled system. In order to realize energy efficient locomotion so that the control inputs can be reduced to a reasonable level and the stability of motions could be improved effectively, it is necessary to carefully plan the target locomotion trajectory of the hopping robot.
A steady continuous motion process of the hopping robots is commonly a periodic motion, which is roughly composed of two harmonic oscillators with different nature frequency, and every motion cycle is composed of four phases, flight phase, touchdown, stance phase, and lift off. During flight and stance phases, the motions of the hopping systems can be controlled, while, in touchdown and lift off phases, the hopping systems can not be controlled due to the short intervals of constraint changes. During stance phases, the nature frequency of the spring leg is so the resonance time period of the linear spring leg is As illustrated in Figure 2, the thick line segment represents the time-amplitude curve of the linear spring leg in stance phase of a motion cycle. Suppose the length of the leg at the instants of touchdown and lift off is the nature length l 0 , the duration of stance phase is T v , the initial phase-angle of the harmonic vibration of the leg is φ l , and the vibration amplitude of the spring leg is l m , then the oscillations of the linear leg in stance phases can be approximately expressed by and the speed of the linear leg is approximately given bẏ wherel(t + vc ) is the compress speed of the linear leg at just after the touchdown. By using the two initial conditions of Equations (64) and (65), the two undetermined parameters l m and φ l can be resolved as From Figure 2, the duration of a stance phase can be calculated as where t φ l = 1 ω l asin mg k l l m . Both in stance phase and flight phase, the vibrations of the hip joint should follow the same harmonic oscillation, but with different initial phase angles. By the dynamics of the hopping robot in flight phases Equation (47), the nature frequency of the hip joint can be calculated as However, the duration of a flight phase is determined by the lift off speed of the robot, and is given as Then, the time period of the hopping motions is Similar to Figure 2, the angular motions of the body and the leg can be analyzed in Figure 3, where the thick line segment denotes the flight phase in a hop cycle. From Figure 3, it is easy to show that the angle trajectory of the leg in the flight phases can be written as and satisfy the boundary conditions where β is the swing amplitude of the leg in the hop cycles and is always given by the control commands. By using the initial conditions Equation (72), the undetermined parameters A and φ θ of Equation (71) can be resolved as for the stance phases t ∈ 0 T v , and for the flight phases t ∈ 0 T f . Since during the flight phases the angular motions of the hopping robots satisfy the law of angular momentum conservation, from Equation (47), the angular speed of the body is determined byφ where we suppose the initial angular momentum of the robot at lift off is zero. The angular position of the body should be and the parameters A and φ θ in Equation (76) are given by Equation (73) for the stance phases, and given by Equation (74) for the flight phases. Note that there already exists a phase difference π in Equations (71) and (76) since the speed directions of the body and the leg are contrary.

Hopping Control of the One-Legged Planar Robot
As shown in Section 3, for a hybrid dynamical system, it is possible to control the system to certain time varying trajectories by partial state feedback and time-sharing patterns. For the one-legged hopping robot system considered in this paper, the time of continuous dynamics modes (stance and flight phases) is very limited, which do not exactly satisfy the condition (A6) of Theorem 3 or (B3) of Theorem 4. Thus, the state jumps caused by the constraint variations of the hybrid systems can quickly destroy the stability of the closed-loop hybrid systems. However, for the one-legged hopping robot system such as that illustrated in Figure 1, we can use the following two level control approaches to stabilize the periodic hopping systems.
During the continuous dynamics modes, linear controllers with only partial state feedback are utilized on the basis of the Theorems 3 and 4. That is, in stance phases, the variables l(t) and ϕ(t) of the robot will be controlled to their target motions, while the rotational motions θ(t) of the leg will be open loop. Thus, the controller in stance phases is where e ϕ = ϕ(t) − ϕ d (t), e l = l(t) − l d (t), t ∈ 0 T v , and k i (i = 1, . . . , 4) are positive constants.
In flight phases, the angle of the leg is only controlled to the given trajectory θ d (t) in closed loop so that the robot lands in right states, while the body is actuated in open loop by the opposite reaction of the leg. In flight phase, the length of the leg does not change, i.e., l(t) = l 0 . Thus, the controller in flight phases is given by where e θ = θ(t) − θ d (t), t ∈ 0 T f ; k 5 and k 6 are also two positive constants. Note that the feedback gains k i (i = 1, . . . , 6) can easily be obtained by the approaches of linear quadratic optimal regulators (LQR). In order to stabilize the one-legged robots to a desired periodic orbit trajectory, an adaptive controller with discrete time feedback is applied to stabilize the leg's attitude in stance phases where k is the current times of touchdown, θ th (k) is the desired angle of the leg at the instant of touchdown, is the change of the leg angle at the instant of lift off, and , r 0 , r 1 and r 3 are adjustable positive constants, and θ 0 is determined by the desired locomotion speed of the robot and is approximately given by θ 0 = π 2 + atan ẋ cm T v 2l 0 . Suppose the physical parameters of the robot model as that given in Table 1, and the initial states of the robot are given by x cm , z cm , θ, ϕ, l,ẋ cm ,ż cm ,θ,φ,l t=0 = 0, 1.0, π 2 − 10π 180 , π 2 + 10π 180 , 0.531, 1.7, 0, 0, 0, 0 , where (x cm , z cm ) t=0 = (0, 1.0) is the initial position of the CM of the robot, (θ, ϕ, l) t=0 = π 2 − 10π 180 , π 2 + 10π 180 , 0.531 is a given initial configuration of the robot, ẋ cm ,ż cm ,θ,φ,l t=0 = (1.7, 0, 0, 0, 0) is the initial speed of the robot. Then, by applying the linear controllers (77), (78) and (79), some numerical simulation results of stabilizing the hopping robot to locomotion with moving speed about 2.0 m/s are illustrated in Figures 4-9. Differential equations are solved by employing the well-known Runge-Kutta algorithm.        Figure 4 shows the trajectories of the variables (θ, ϕ, l) in the last three cycles of fifty hops in total, Figure 5 illustrates the corresponding control inputs, and Figure 6 plots the animations of the robot in the last three hops. The state jumps measure Equation (60) of the controlled robot system is illustrated in Figure 7, which shows that the state jumps at touchdown phase finally converge to zero. In other words, the hopping robot is stabilized to a periodic orbit with zero impulse at touchdown phases, so that the energy efficiency of the hopping system is improved and the maximum of the inputs is reduced to a reasonable level (referring to Figure 5); otherwise, the maximum of the inputs will be very large under closed-loop mode. However, the nonzero inputs of the robot in stable hops are caused by the approximate target trajectory Equation (64) since the spring leg in stance phase is actually an oscillating system with variant stiffness. Referring to Equation (52), the actual equivalent stiffness of the linear spring leg should bek l = k l − mθ(t) 2 . In addition, as examples, Figures 8 and 9 also respectively show that the state variables l(t),l(t) and θ(t),θ(t) are stabilized to a period one orbit.

Conclusions
From the viewpoint of Lagrange-d'Alembert equations with additional constraints, this paper investigates the modeling approaches of a class of VCS with instantaneously variant constraints, and the switching control techniques of stabilizing the VCS to given periodic orbits. We show that under certain conditions the VCS are essentially a class of HDS, and there possibly exist periodic orbits with zero impacts, so that the HDS can be stabilized to period-one orbits by linear controllers with partial state feedback, even though the HDS are generally underactuated nonholonomic systems in time continuous modes. These new discoveries point out that: (a) in an HDS, there generally exist time-varying trajectories without impacts, and, (b) if the zero impact trajectories are physically meaningful, then the controlled HDS can be stabilized to the zero impact trajectories by simple switched linear controllers, and (c) the zero impact trajectories provide high energy efficient locomotion modes since the energy loss caused by impacts is avoided.
In addition, by a one-legged planar hopping robot, the main results of modeling and control approaches for the HDS are discussed in detail, and the numerical simulation results are in good agreement with the conclusions of this paper. In order to use the proposed control approach in practice, improving the robust stability of the closed-loop system and applying the pulse control [30] should be further considered.

Conflicts of Interest:
The authors declare no conflict of interest.

VCS
Variant constraint systems HDS Hybrid dynamical systems DOF Degrees of freedom CM Center of mass