A Tutorial for the Analysis of the Piecewise-Smooth Dynamics of a Constrained Multibody Model of Vertical Hopping

Contradictory demands are present in the dynamic modeling and analysis of legged locomotion: on the one hand, the high degrees-of-freedom (DoF) descriptive models are geometrically accurate, but the analysis of self-stability and motion pattern generation is extremely challenging; on the other hand, low DoF models of locomotion are thoroughly analyzed in the literature; however, these models do not describe the geometry accurately. We contribute by narrowing the gap between the two modeling approaches. Our goal is to develop a dynamic analysis methodology for the study of self-stable controlled multibody models of legged locomotion. An efficient way of modeling multibody systems is to use geometric constraints among the rigid bodies. It is especially effective when closed kinematic loops are present, such as in the case of walking models, when both legs are in contact with the ground. The mathematical representation of such constrained systems is the differential algebraic equation (DAE). We focus on the mathematical analysis methods of piecewise-smooth dynamic systems and we present their application for constrained multibody models of self-stable locomotion represented by DAE. Our numerical approach is demonstrated on a linear model of hopping and compared with analytically obtained reference results.

Alternative modeling approaches do not necessitate the experimental results, but develop predictive models, which imitate the dynamics of locomotion. For instance, references [8,9] introduce models for the investigation of foot impact intensity at forefoot, mid-foot and rear-foot ground collision. In [10], the geometric non-linearity of a static segmented leg model is studied. The SLIP model in [11] has been developed for the understanding of the essential dynamics of running and hopping locomotion. The further detailed variation of the SLIP model is presented in [12]: the segmented leg provides a geometry very similar to the human leg; however, the inertia of leg segments is neglected; consequently, the ground-foot collision and the impact-induced energy absorption are not investigated. Reference [13] introduces a novel approach for humanoid gait generation. The control of foot placement and human balancing problems are discussed in [14], where a predictive model is applied.
Our research contributes to the development of mathematical models of legged locomotion. Since the mathematical description of the above-mentioned and the similar predictive models are quite challenging, these models develop slowly, and the integration of the features leads to highly intricate mathematical description. The long-term goal of our research is to reach such integrated models of hopping and running that can analyze (i) the dynamic stability issues, the neural processes such as balancing [14][15][16][17], the motion control and the motion pattern generation [18]; (ii) geometric non-linearities [10,12,19,20]; (iii) ground-foot collision intensity and impact-induced energy absorption [9,21,22]; (iv) distinct topology and actuation in the flight and ground phases [23].
In paper [24], we introduced a simple 2 DoF model of hopping that demonstrated how the active control provides stable hopping with uniform hopping height despite of the impact-induced energy absorption. The model consists of two mass blocks one above the other and connected to each other by means of an active spring-damper. In [24], we demonstrated a general dynamic analysis method, which is suitable for more complex models, such as the 6 DoF model in [20]. The analysis of the three-segmented hopping model in [20] is based on the dynamic analysis method explained in [24] and demonstrates the applicability for higher DoF systems.
The present paper introduces and overview a dynamic analysis methodology which is suitable for higher complexity multibody models. The mathematical representation of the models in [20,24] are based on the classical general coordinates of the minimum set. However, in the case of high number of DoF or closed kinematic loops, the application of non-minimum set (also referred as redundant set) of descriptor coordinates is more efficient. The application of non-minimum coordinates in multibody dynamics is well known from [25][26][27][28][29][30]. In case on non-minimal coordinates, the dynamic equation of motion is subjected to algebraic constraint equations that assure the geometric connection of the redundant coordinates. The constrained representation of the equations of motion suits well the possibly intricate combination of ground-foot contacts and makes possible to analyze models when more than one leg can touch the ground. This technique leads to the situation, when the equation of motion has DAE form. In the present work, we combine the stability analysis method of the periodic solution of piecewise-smooth dynamic systems that are described by ordinary differential equations (ODE) [21,[31][32][33] with the efficient DAE formulation of multibody systems [25,27,34] and legged locomotion systems [28,30]. The fundamental theories and ideas of the analysis of the legged locomotion with similar approaches are detailed in [35]. The applied dynamic analysis technique is specialization of the method in [36] for hopping and running multibody models.

Organization of the Paper
Our models of vertical hopping with coordinates of minimum set and redundant set are detailed in Section 2. The models consist of the equation of motion of the mechanical system and the equations of the control actions. Section 3 details the dynamic analysis method, which is specifically tailored for non-linear hopping and running systems with redundant coordinates and two continuous phases. The numerical results together with their analytical validation are presented in Section 4. The analytical calculation exploits the fact that the model is linear; however, it is not a requirement for the presented numerical dynamic analysis method. The conclusions are summarized in Section 5.

A Constrained Model of Vertical Hopping
The DAE-based dynamic analysis methodology of piecewise-smooth hopping and running systems is presented by the case study of an excessively simplified model of hopping. Hence, we can rather focus on the analysis methodology.
For the sake of simplicity, a constrained linear model of hopping is introduced in this section. However, non-linearities may emerge in locomotion models due to geometric non-linearities, such as in case the model in [37]; spring characteristics due to the muscle-tendon dynamics [10]; or non-linear behavior of the neural model [11]. Although these non-linearities do not appear in our basic hopping model, the dynamic analysis method is not limited to linear systems.

Mechanical Model
Our simplest hopping models consist of two vertically moving particles connected by an active spring-damper. The model is placed in gravitational field above smooth, rigid, and horizontal ground. The constrained model, which is shown in the right panel of Figure 1, has been developed based on the unconstrained hopping model that was published in [24] and is shown in Figure 1 left. The unconstrained model in [24] and our new model match each other in terms of physical parameters and DoF. The difference is the mathematical representation: the upper and lower body has been split into two parts and constrained to each other in the new model. The constrained model on the right panel of Figure 1 is used as a case example to illustrate our analysis methodology; however, the presented methodology is suitable for multibody models with high complexity. Similarly to the unconstrained model, the physical parameters are the total mass m and the ratio µ ∈ [0, 1] of the upper and lower mass. The upper m U and lower masses m L are further separated by means of parameters µ U ∈ [0, 1] and µ L ∈ [0, 1], with which where are the mass of the upper and lower blocks in the unconstrained model. Further parameters of the model are the stiffness k and the damping parameters d F and d G , which are valid in flight and ground phases, respectively. For the better interpretation of the results (in Section 4) and in order to reduce the number of parameters, we introduce α = √ k/m, D F = d F /(2αm) and D G = d G /(2αm).
The system is described by the non-minimum set of coordinates q = [z 1 (t), z 2 (t), z 3 (t), z 4 (t)] T of the number of n = 4, and the following geometric constraints ensure the connection of the descriptor coordinates: The rigid body constraints in γ b are permanently valid in both the flight and the ground phases. The number of permanent rigid body constraints is h = 2 in our case example. In the case of this hopping model, γ b is linear. However, it is not a requirement for the general formalism detailed in Section 3.2.
In the specific case of our hopping model, the number of constraints that represent ground contact is l = 1, and the corresponding constraint vector is Again, in the case of our model, γ g is linear too, which is not a requirement for the general formalism detailed in Section 3.2.
To construct the equation of motion, we derive the mass matrix H and the vector C of spring, damping and gravitational forces, which are where L 0 is the unstretched spring length. The general force vector C in (10) and the mass matrix in H (9) are linear and constant, respectively. However, C and H can be non-linear functions of the state variables in the general formalism (Section 3.2).

Active Spring-Damper for Ensuring Periodic Motion
The flight (airborne) phase and ground (support or stance) phases, which periodically follow each other, are clearly separated in our models. In such systems, some segments of the limbs periodically come into contact with the ground; therefore, the motion of these body segments are constrained during the whole duration of the ground contact. The ground collision at the end of the flight phase leads to the loss of a certain amount of kinetic energy [9,11,21,24]. This energy amount is referred as constraint motion space kinetic energy (CMSKE) in [22]. As the muscles or actuators provide mechanical power, the kinetic energy accumulates in the rest of the body during the ground phase. In the case of periodic motion, the mechanical energy reaches its original level at the end of the ground phase. The ground contact then terminates, and the flight phase begins again.
The unconstrained model, which is depicted in the left panel of Figure 1, is developed for the illustration and dynamic analysis of the above explained sequence. The ground collision of m L terminates the flight phase and the total amount of the kinetic energy of m L is absorbed, since the collision is modeled to be fully inelastic. The active spring-damper plays the role of the muscles: it accelerates m U in upward vertical direction and increases its kinetic energy. If the height of the new hop is the same as it was in the previous one, then the motion is periodic.
We can conclude that hopping systems with collisions are non-conservative without active control. This issue and the modeling of energy supply are thoroughly surveyed in [11]. A specific concept for the energy refilling of impacting hopping systems is presented in the case of a general 1 DoF non-linear model of hopping presented in [38]. Since the damper and the spring is linear in our model, it is more specific. Having linear dynamics allows us to carry out analytical calculations, with which the numerically obtained results are validated (see Section 4.1). Besides, the presented numerical analysis is not limited to linear systems.
The active spring-damper system imitates the behavior of limb muscles in the sense that muscles supply energy lost to dissipation and foot impacts according to [11,38]. The spring-damper plays different role in the flight and the ground phases. The damper d has a certain positive value d F in the flight phase to suppress undesired oscillations. However, it has a properly chosen negative value d G in ground phase to regenerate mechanical energy. The negative damping could be achieved by means of an active actuator in a robotic realization of the model. The stiffness k is constant regardless the variation of the flight and the ground phases. The operation of the active spring-damper system is the same in the constrained model.

Switchings between the Flight and the Ground Phases
As we explained in Section 2.2, the flight phase terminates when the mass block in the lowest position reaches the ground. The event of flight to ground (F2G) transition is mathematically expressed by means of the indicator function h F2G ∈ R: The indicator function h F2G crosses zero from positive to negative direction at touchdown. The F2G transition not only cause the switching of the governing equation of the continuous dynamics, but involves an impact too. The impact is modeled as a mapping (i.e., resetting) of the state variables: the velocity of the lower mass blocks becomes zero and other state variables remain the same. The mapping is defined as The ground to flight (G2F) transition occurs, when the ground-foot force becomes a pulling force. Mathematically, the G2F transition occurs, when indicator function crosses zero from negative to positive direction. The Lagrange-multiplier λ g in h G2F ∈ R is the ground contact force acting on particle m 4 . Negative and positive sign of λ g refers to pushing and pulling force, respectively. The mapping after ground-foot release is identity, since the state variables does not change. This type of transition characterizes the Filippov systems [31,39]. The mapping is defined by

Dynamic Analysis
The analysis methods of the periodic solutions of hybrid systems, which are described by minimum coordinates, are detailed in [21,31,32,39]. In particular, the analysis and control of the local stability of periodic orbits of hybrid systems is detailed in [32]. A constrained formalism for the analysis of legged locomotion systems is developed in [28,30]. The general analysis of piecewise constrained systems is detailed in [36]. The formalism is particularly efficient when intricate combination of ground-limb contact emerges. Our contribution is in the application of the stability analysis method in multibody models of hopping and running systems in constrained mathematical representation. Furthermore, the method is demonstrated on an illustrative case example of a simple hopping system.

Assumptions and Limitations
We focus on hopping or running motion only, when the limb(s) that periodically touches the ground has a larger than zero elevation in the flight phase. Accordingly, motion types with continuous ground contact are out of scope in this work. We assume that (a) the flight phase and the ground phase periodically switch each other. This assumption helps to avoid the consideration of intricate combinations of the sequence of various dynamics, which may happen e.g., in the case of more than one limb.
Furthermore, we assume that (b) the ground-foot collision is instantaneous, which leads to infinitely large instantaneous forces over infinitesimal time duration so that the net impulse due to the impact force is finite [8,9,22]. Additionally, (c) completely inelastic collision is assumed, hence there is no rebound [8,22]. Assumption (d) is that there is no slip (more details in Section 3.2). Assumptions (b), (c) and (d) allows us to consider the ground-foot contact as an instantly arising geometric constraint such as in [1,9,22,[40][41][42].

Continuous Dynamics in the Flight and in the Ground Phases
The general form of the equation of motion is a DAE, since the dynamic Equations (15) are augmented with the algebraic equations that represent geometric constraints (16): where the vector q contains the differential variables and the vector λ contains the algebraic variables. The mass matrix of the unconstrained system is H ∈ R n×n , the vector C ∈ R n represents the velocity dependent inertial forces, the forces exerted by springs and dampers, gravitational forces and any active forces in general. The matrix is the gradient of the geometric constraint vector γ, and λ is the vector of the corresponding Lagrange-multipliers. The equation of motion (15) together with the geometric constraint Equation (16) are reorganized by means of DAE index reduction, which leads to the following general form used in [25,27,28,30]: The form (18) makes possible to express the accelerationq ∈ R n as the function of the state variables, which are collected in x = [q,q] T , where x ∈ R 2n is the state vector. Finally, the continuous dynamics is governed byẋ where f F ∈ R 2n and f G ∈ R 2n are smooth vector fields corresponding to the flight and the ground phase dynamics, respectively.
The vector field f F is obtained by using the permanent rigid body constraints γ b ∈ R h only (for our case example, see (7)). Hence, the entire constraint vector γ F ∈ R h is the following in the flight phase: The ground-foot contact is modeled by geometric constraints γ g ∈ R l . In our specific case example (8) stands. In the ground phase, the overall constraint vector γ G ∈ R h+l , which is incorporated in f G , reads: Considering a convenient running kinematics, assumption d in Section 3.1 is usually satisfied in practice, since depending on the surface quality, normally there is no significant slip of the foot. Although the slip-free condition is checked in such models, when not only vertical motion is considered, e.g., [19,20]. The ground-foot constraint is separated into two parts: γ t and γ n respectively constrains the motion in the tangential and in the normal direction relative to the ground surface. The ground-foot constraint reads: in the frame of the tangential and the normal directions. To satisfy the non-slip condition, the Coulomb friction coefficient must be larger than the critical value [43]: where γ q,t and γ q,n are the Jacobians of constraints γ t and γ n respectively.

Piecewise-Smooth Periodic Orbits and Impacts in Hopping and Running
The combined solution φ(x 0 , t) of (19) and (20) go through the smooth vector fields f F and f G , which are separated from each other by switching surfaces Σ F2G and Σ G2F as it is shown in Figure 2. Flight to ground (F2G) transition occurs at t = t F2G , when the switching surface Σ F2G , which is defined by the scalar indicator function h F2G (x), is reached from positive to negative direction. The indicator function h F2G (x) expresses the relative position of the ground and the end of the limb. After the impact, the solution remains on the Σ F2G switching surface itself in the entire ground phase. Ground to flight (G2F) transition occurs at t = t G2F , when set Σ G2F , which is defined by h G2F (x) and intersects Σ F2G , is crossed from negative to positive direction. Function h G2F (x) expresses the normal component of the ground contact force. We remark, that the entire period ends at time instant t p := t G2F .
According to the PhD thesis by [31], piecewise-smooth systems can be classified according to the degree of discontinuity. In the case of the G2F transition, the solution goes through the switching surface Σ G2F without discontinuity (i.e., without mapping or reset of the state variables). Hence, the mapping is an identity: However, the vector field changes abruptly and therefore the solution is non-smooth continuous, such as in the case of Filippov-type systems.
On the contrary, F2G transition involves an inelastic impact, which cause the discontinuity of the solution: the velocity of the contact point(s) becomes zero abruptly when the solution reaches the switching surface Σ F2G , i.e., when the limb reaches the ground. The jump function g F2G (x) maps the state into a specified location in the state space (see [32]). In general, the positions do not change, and the velocities are projected by P a ; therefore, the jump function yields where P a projects into the admissible direction of the set of constraints γ G defined in (22). References [25,30,34] give the formula for the calculation of the projection matrix where γ G,q is the Jacobian of γ G . Piecewise-smooth systems with impacting discontinuities are referred as hybrid systems in [21]. The projection is applied for the calculation of the constrained motion space kinetic energy (CMSKE) too, which vanishes at each ground-foot impact. The admissible motion space kinetic energy T a and CMSKE denoted by T c are calculated with the pre-impact general velocity vectorq as [43] T a = 1 2q T P T a HP aq , with P c = I−P a . Papers [22,43] showed that foot strike intensity can be characterized by the CMSKE, which depends on the pre-impact configuration and velocity and the effective mass matrix H e = P T c HP c . The effective mass concept for foot impact was introduced first by [44] for a one DoF model of legged locomotion. According to references [22,43], CMSKE is directly proportional to the impulse of the contact reaction force and to the peak reaction force; therefore, CMSKE is used as an indicator of the foot impact intensity in our study.

Numerical Identification of the Periodic Orbits
The periodic orbits, sketched in Figure 2, are found by means of the shooting method: the initial conditions x 0 = x(0) at the beginning of the period are varied, until the solution at the end of the period coincides with the initial values.
The starting time instant t = 0 of the period is chosen to be the beginning of the flight phase. The advantage is that a subset of the state variables are known here, since we initiate the motion from switching surface Σ G2F ; i.e., the motion starts at the time instant when the ground-foot contact releases. The variable subset of initial values x 0 are collected in subvector u 0 (x 0 ) ∈ R u . The periodic orbit is found by the shooting method as the root x * 0 of the problem where u e is the corresponding subvector at the end point of the solution trajectory. We note that subvector u 0 has lower dimensions than x 0 ; therefore, the number of unknowns is reduced, which causes considerably decrement in the computation time. In the specific case of our hopping model, the subvector of unknown state variables yields: with u = 6, while z 4 = 0 andż 4 = 0 are set as initial conditions at t = 0.

Stability Analysis
We use two alternative approaches for the determination of the stability of the periodic orbit. A computationally cheap option is to gain information about the stability based on the approximation of the Jacobian M(u 0 ) ∈ R u×u of the solution with respect to the initial values at the end of the period: The computational efficiency of this approach is induced by the fact that the numerical approximation of the Jacobian J(u 0 ) = ∂F(u 0 )/∂u 0 is necessarily calculated in the shooting method. Based on (30) and (32), the approximation of M(u 0 ) is gained as M(u 0 ) = J(u 0 )+I, (33) where I is the identity. As a stability criterion we can say that the periodic orbit is stable, if a perturbed initial condition u 0 is mapped "closer" to the reference (periodic) path. Thus, ∀ |µ i | < 1 is the criterion of asymptotic stability, where µ i are the eigenvalues of M. The accuracy of M might not be satisfying in the case of more complex systems due to the discontinuities. Instead, a computationally more demanding, but at the same time more accurate computation approach of the Jacobian of the solution is suggested by [21,31,39,45]. The coupled ODEs (34) and (35) of size 2n+(2n) 2 is solved along continuous phases: where I is 2n×2n identity matrix and f x is the gradient of the smooth vector field f. Equation (35) is referred as first variational equation in the literature. Equations (34) and (35) are applied with substitutions f = f F and f = f G in flight and ground phases, respectively. The Jacobians Φ F (t) and Φ G (t) (referred as fundamental solution matrix in [31,39]) of the trajectories φ F (x 0 , t) and , t) are obtained in the flight and ground phase, respectively. The JacobianΦ(t) of the composite trajectory φ(x 0 , t) along the entire period contains the Jacobian S of the discontinuity mapping too, which is referred as saltation (or jump) matrix in [21,31]. The monodromy matrix, i.e., the composite flow Jacobian at the end of the period yields The monodromy matrixΦ(t G2F ) can be reliably applied for judging asymptotic stability of periodic orbits, by using criterion ∀ |ν i | < 1, where ν i are the eigenvalues ofΦ(t G2F ) (referred as Floquet multipliers in the literature). Accordingly to reference [45,46], there is always one eigenvalue, which equals to 1, since f(x(0)) = f(x(t p )) is true for the vector field in (34).
We remark thatΦ(t) can be used as a Jacobian of the trajectory where φ(x 0 , t) is the perturbed discontinuous trajectory, φ(x ref 0 , t) is the reference trajectory andΦ(t) is the state flow Jacobian having discontinuities too.
The saltation matrices in (36) are calculated based on the literature. First let us consider the G2F transition (ground-foot release), where there is no discontinuity mapping of the solution, as it yields from (25). The saltation is given by the formula published in [31]: where ) and h G2F,x = ∂h G2F /∂x is the gradient of the indicator function evaluated at x(t G2F ). The discontinuity mapping g F2G from (26) is taken into account in the F2G transition (ground-foot collision). The corresponding formula, which is a generalization of (38), is obtained from [21,32]: F2G,x = ∂h F2G /∂x and g − F2G,x = ∂g F2G /∂x are the gradients of the indicator and jump functions respectively evaluated at x(t − F2G ). Here, the time instants t − F2G and t + F2G are right before and right after the foot impact.

Results and Validation by Means of Semi-Analytic Calculations
This section presents the results obtained for the models depicted in Figure 1. To obtain a reference, with which the numerical results are compared, a semi-analytic calculation is carried out for the unconstrained model. Then the results for the constrained model are presented: first, the periodic orbit is shown at a certain parameter set; and secondly, the effect of physical parameters are studied.

Semi-Analytic Calculations
The dynamics of the unconstrained model (see left panel in Figure 1) is analyzed to validate the numerical results. First, the closed form solution for the flight and ground phase is obtained and then the periodic orbits are generated by directly solving the non-linear algebraic equations that express the boundary conditions.
For the sake of simplicity, L 0 = 0 is applied in the analytical calculations. This simplification causes the shifting only of the z U (t) trajectory.

Closed form Solution in the Flight Phase
The analytical solution for the linear equations of motion is obtained using the modal transformation of the general coordinates according to [47,48]. The modal transformation requires the following linear equation of motion for position coordinates z = [z U , z L ] T : with the constant mass matrix H, damping matrix D, stiffness matrix K and the constant vector g of external forces. The system matrices are the following: By solving the characteristic equation det(−ω 2 n H+K) = 0 for the eigenfrequencies ω n,i and by solving the linear equation (−ω 2 n,i H+K)A i = 0 for the modal shapes A i , we obtain ω n,1 = 0 ; The modal transformation matrix T = [A 1 , A 2 ] gives the vector of modal coordinates as ζ = T −1 z. The modal coordinates are The modal coordinate ζ 1 corresponds to the center of mass position of the entire model and ζ 1 corresponds to the position difference of the mass blocks. The solution for these coordinates is a simple parabolic function caused by falling and an exponentially decaying harmonic oscillation respectively, which is shown by the transformed equations of motion for ζ. The transformed equation of motioñ Hζ+Dζ+Kζ+g = 0 withH = T T HT,D = T T DT,K = T T KT andg = T T g yields The general solutions for the differential Equations (45) and (46) are with unknown constants b 1 and b 2 , with undamped natural frequency ω n = ω n,2 , damping ratio ξ = d F (m U +m L )/(2m U m L ω n ) and damped natural frequency ω d = 1−ξ 2 ω n . After the inverse modal transformation z = Tζ, the closed form solution for the original position coordinates z U = ζ 1 +m L ζ 2 and z L = ζ 1 −m U ζ 2 are: To obtain simple form of the boundary conditions, the integration constants b 1 , b 2 , b 3 and b 4 in (49) and (50) are expressed in terms of the initial conditions z U,0 = z U (0), z L,0 = z L (0),ż U,0 =ż U (0) andż L,0 =ż L (0).

Closed form Solution in the Ground Phase
The mass block m L is fixed to the ground during the entire grounded phase; therefore z L (t) = 0 stands. The governing equation of the 1 DoF ground phase dynamics and its closed form solution reads with undamped natural frequencyω n = √ k/m U , damping ratioξ = d G /(2m Uωn ), damped natural frequencyω d = 1−ξ 2ω n and unknown integration constantsb 1 andb 2 .
For the sake of simple expression of the boundary conditions, the shifted time axis τ is introduced as where t F is the time duration of the flight phase. The ground phase time duration and the entire time period are τ G and t p = t F +τ G respectively.

Periodic Solutions
The closed form solutions for the continuous phase dynamics obtained in Sections 4.1.1 and 4.1.2 makes possible to get rid of the computational costs and inaccuracies of the numerical time integration. The initial conditions that correspond to the periodic paths are obtained as the solution of the following algebraic equation system: where (54) and (55) assures that the motion of the mass block m L in the beginning of the flight phase is initiated from the ground with zero velocity at t = 0. Equations (54) and (55) gives the trivial solutions z L,0 = 0 andż L,0 = 0; therefore, the number of equations is reduced from 8 to 6. Conditions (56) and (57) guarantees the continuity of the position and velocity of mass block m U in the F2G transition (ground collision). The indicator function h F2G (x) , which is satisfied at the F2G transition, is expressed in (58). The continuity of the position and velocity of the mass block m U is assured by (59) and (60). The indicator function h G2F (x) is expressed by (61). The contact force λ(τ) is a simple expression of z U andż U , since z L =ż L = 0 in ground phase [24]: The non-linear algebraic Equations (56)-(61) are solved for the unknown variables by means of an iterative solver (such as the Newton-Raphson method) at the different values of the bifurcation parameter t F . The reason for choosing t F as bifurcation parameter is that every other parameter and indicator are strictly monotonic function of t F ; therefore, the numerical issues in fold points are avoided. Since the non-linear Equations (56)-(61) solved numerically, the results are semi-analytic.

Stability Analysis
In the present analytical approach, the stability is judged based on the Jacobian M of the solution with respect to the proper initial conditions. The Jacobian is calculated directly. The initial conditions of the reference trajectory φ(z 0 , t), which corresponds to the periodic orbit, were perturbed by a small number ∆ = 10 −4 . Then two trajectories φ(z 0,∆1 , t) and φ(z 0,∆2 , t) were obtained in the form of (49), (50) and (52) with the perturbed initial conditions respectively. Note, that these initial conditions are off the G2F switching surface, from which the reference trajectory starts. As Figure 3 shows, the flight phase trajectories propagate separately from the reference trajectory. The F2G transition occurs in different points of the Σ F2G switching surface and in different time instants compared to the reference trajectory. Similarly, the G2F transition, which is the end of the entire period, occurs in separate points of Σ G2F in separate time instants. The Jacobian M ∈ R 2×2 is obtained as where φ U is the trajectory of the z U (t) state variable.

Illustration of Periodic Solutions
This section shows the periodic solution for a case example with a certain parameter set, which is the following:  Figure 3 shows a typical periodic orbit φ(z 0 , t) in the z 1 , z 2 andż 1 subspace. Also, the perturbed trajectories φ(z 0,∆1 , t) and φ(z 0,∆2 , t) are illustrated. For the sake of better illustration of the idea of the stability analysis, the initial conditions are perturbed by a large number instead of the small number ∆ (see Section 4.1.4). The perturbations of initial conditions (65) and (66) are therefore visible. The quantities φ(z 0,∆1 , t p )−z 0 and φ(z 0,∆2 , t p )−z 0 , from which the Jacobian matrix (67) is constructed, are also visible here.
The stability is determined in two alternative ways, as we explained in Section 3.5. The monodromy matrixΦ (see 36) was calculated by means of the variational Equations (34) and (35). Additionally, the approximation of the Floquet multipliers was calculated with the help of the Jacobian M (see (32) and (67) in the specified parameter set. The system is stable with the current parameter set, since the largest relevant eigenvalue is µ 1 = ν 2 = 0.4714 < 1.
To illustrate the discontinuities in the system, the sections of the phase-space are shown in Figure 4 separately for the coordinates z U and z L . The jump defined by the function, which appears in (12) and (26), is illustrated in the right panel. The switching surfaces Σ F2G defined by (11) and Σ G2F defined by (13) are visible in the right and left panel of Figure 4 respectively. Figure 3. The switching surfaces, the periodic orbit and the perturbed trajectories. The dots indicate the intersection with switching surfaces Σ G2F and Σ F2G . Note: state variableż 2 has a jump after reaching Σ F2G , which cannot be seen explicitly here. Figure 5 shows the time history of the mechanical energy. The constrained motion space kinetic energy T c and the admissible motion space kinetic energy T a (CMSKE and AMSKE) can clearly be seen here. The kinetic energy of the lower mass block is constant zero during the entire ground phase because of the inelastic ground-foot collision. Although CMSKE is absorbed in each F2G transition, the energy level at the end of the period reaches the value that corresponds to the beginning of the period. The extra energy is pumped to the system by the negative damper during the ground phase.
The time histories for all state variables were in correspondence in the semi-analytical and in the numerical calculation.

Overview of Effect of Parameters
To gain an overview of the effect of physical and control parameters, the portion µ of the masses and the negative damping value D G were varied by means of a numerical continuation. Parameter µ was set to the discrete values 0.8 and 0.95. For each µ values, D G was continuously decreased from zero until the periodic orbits disappear. Compared to [24], it is a new result about the hopping model that the unstable periodic orbits were found, not only the stable ones.
Plotting the flight phase time period t F , the ground phase time period t G and the entire time period t p in Figure 6 provides a general look of the dynamic behavior of the system. Above D G = 0 values there is no periodic solution, since both dampers continuously absorb energy, which leads to static equilibrium. In the D G = 0 case, the time period t G of the ground phase is exactly equal to the time period t p0 = 2π √ m u /k of the 1 DoF oscillatory system composed by spring k, damper d G = 0 and mass m u . For D G > 0 cases, the ground phase time period t G becomes infinity, since the lower mass block never takes off from the ground. When increasing D G in the negative direction, one or more stable periodic orbits can be found with different time periods. Finally, the flight phase time period goes to infinity. This is the point, when the Floquet multipliers shown in Figure 7 would cross the 1 value.   Figure 8 shows the apex height of the z U and z L trajectories. The hopping height is in close relation with the flight phase time duration t F , which can be seen by comparing the graphs in Figures 6 and 8. In the extreme values of the ground phase damping parameter D G , the hopping height is such small or such large that are practically irrelevant. In the middle range, the hopping height is realistic for humans, according to [6,7]. The ratio of the CMSKE and the AMSKE (absorbed and preserved amount of kinetic energy) to the total kinetic energy T are plotted in Figure 9. One can notice that the AMSKE value denoted by T a converges to the µ value as D G decreases. We can conclude that it is possible to find optimal parameter set when the ratio of CMSKE is the smallest. Since CMSKE in in relation to the impact intensity and is indirectly in relation to the risk of injuries, this may be a cost function in running and walking systems. The Floquet multiplier, which determines the stability of the periodic orbits, is plotted in Figure 7. The results show that the stability highly depends on both physical (µ) and control (D G ) parameters. At D G = 0 the Floquet multiplier is exactly 1, since in that case a 1 DoF undamped oscillation occurs with constant amplitude. By decreasing the numerical value of D G , the Floquet multiplier reaches 1 and the time period of the orbit becomes infinitely long (see Figure 6). Comparing Figures 7 and 8, one can notice that stable periodic orbits are separated from each other by unstable periodic orbits. The combination of apex height, stability and impact intensity may be applied as a combined cost function in certain applications in the field of legged robotic systems.
Our dynamic analysis methodology presented in Section 3 makes possible to obtain similar continuation result for complex multibody models with high DoF and geometric constraints.

Conclusions and Future Work
We further developed a minimally complex model of hopping resulting in a constrained model, on which our dynamic analysis methodology was demonstrated. In general, the constrained models are described by a non-minimum set of coordinates and the equations of motion are augmented with geometric constraint equations. This modeling approach, which is typical in multibody dynamics, allows us to build high DoF models effectively. The dynamic analysis methodology of piecewise-smooth systems was generalized for such constrained multibody systems in this study.
The presented dynamic analysis methodology focuses on the stability investigation of periodic orbits of hopping and running models that consider control issues, geometric non-linearities, inelastic ground-foot collision and energetic costs in relation to the collisions. The controller in our models generates the motion pattern without any prescribed kinematics. Instead, such sets of the physical and the control parameters are identified that assure stable periodic motion.
The limitation of the model and the presented dynamic analysis is that the motion is strictly divided into two phases: flight and ground. Furthermore, we consider only single support in the ground phase. The dynamic analysis methodology will be further developed to handle multiple ground support.
We demonstrated our dynamic analysis method on a vertically hopping constrained model. The results were compared with semi-analytically obtained reference results of the corresponding unconstrained model: the periodic orbits, the parameter set that guarantees stable operation, the Floquet multipliers and energy loss induced by the ground collision were in total correspondence.
In future work, we will apply the present dynamic analysis for higher DoF models of hopping and running. We will particularly focus on ground-foot collision, since foot collision is the main source of injuries when walking and running. Furthermore, a considerable amount of energy-demand of locomotion is related to ground-foot collision.
Author Contributions: R.R.Z. wrote the paper and created the graphical content; B.B. accomplished the analytical calculations presented in Section 4; L.B. contributed in the literature study; furthermore, as a corresponding author, L.B. managed the submission of the paper; A.Z. developed the models and the code for the numerical dynamic analysis presented in Sections 2 and 3 respectively.