Nonholonomic dynamics of the Twistcar vehicle: asymptotic analysis and hybrid dynamics of frictional skidding

The Twistcar vehicle is a classic example of a nonholonomic dynamical system. The vehicle model consists of two rigid links connected by an actuated rotary joint and supported by wheeled axles, where nonholonomic constraints are assumed to impose no skidding of the wheels. Recent experimental measurements conducted with a robotic Twistcar prototype have shown disagreements with previous theoretical analyses. In particular, significant skidding has been observed, in addition to discrepancies with respect to theoretical predictions of divergence in oscillations of the vehicle’s speed and orientation, as well as direction reversal depending on the vehicle’s structure. The goal of our research is to resolve this disagreement by generalizing the theoretical analysis. First, we extend previous asymptotic analysis by incorporating the effects of links’ inertia and oscillation amplitude of the input angle on the direction of net motion. Next, we formulate the vehicle’s hybrid dynamics under frictional bounds and skid-state transitions. Using numerical analysis, we obtain optimal values for the vehicle’s mean speed and energetic cost-of-transport as a function of the input frequency. Our results improve the agreement between theory and experiments and suggest directions for further experimental investigation.


Scientific background
Underactuated locomotion systems are mobile robots whose number of degrees-of-freedom is larger than the number of controlled inputs. Thus, they typically contain passive degrees-of-freedom. The dynamics of such systems is often governed by nonholonomic constraints, which are non-integrable and locally restrict the system's feasible directions of generalized velocities [1,2]. A classical example is Chaplygin's sleigh [3,4] which is a planar rigid body with a blade which is restricted to slide along a body-fixed direction. The same principle holds for vehicles with unactuated wheeled axles which are restricted not to skid (i.e., slip in lateral direction).
In many of these nonholonomic systems, one assumes "kinematic inputs" directly prescribing m actuated joint angles or velocities, which indirectly affect the motion of all N degrees-of-freedom through interaction with the constraints. If the number of constraints n equals N − m, then the system is termed "purely kinematic." Examples of such systems are the Dubins' car [5], three-wheel kinematic snake [6] and truck-trailer(s) chains [7,8]. On the other hand, if n < N − m, then the system's motion involves dynamics, and is governed by interaction between the nonholonomic constraints and balance of generalized momentum variables. Examples of such systems are the snakeboard [9,10], roller racer [11,12], bicycle-like car models [13,14], and several actuated variants of Chaplygin's sleigh [15][16][17]. There have been many works on dynamics and control of nonholonomic locomotion systems, exploiting their geometric structure for analyzing periodic control inputs [18,19].
Another closely related type of robotic locomotion systems are articulated swimmers, where velocity relations that are analogous to nonholonomic constraints are induced by the fluid-swimmer interaction. For micro-swimmers, where inertial effects are negligible, these relations originate from imposing quasistatic equilibrium [20,21], whereas they come from addedmass effect and momentum conservation for inertiadominated swimming in ideal fluid [22,23]. Several works have analyzed such systems utilizing methods of differential geometry [24], optimal control [25][26][27] and perturbation expansion [28][29][30].
For wheeled vehicles moving on low-friction terrain, the no-skid nonholonomic constraints are not always realistic. For example, the kinematic three-wheel snake skids upon reaching singularities at symmetric configurations, as observed both experimentally [6] and theoretically [31]. Several relaxations of the no-skid constraints have been proposed and analyzed. The "skidangle" model [32][33][34][35] introduces skid velocities with non-dissipative constraint forces. In contrast, assuming skidding under viscous frictional forces [36,37] also involves energy dissipation. While those two models involve continuously varying friction forces and skid velocities, Coulomb's dry friction model may induce non-smooth behavior of dynamic skid-state transitions [31,38,39] that also dissipate energy.

Twistcar-Review of recent relevant works
We now introduce the Twistcar toy and its simple theoretical model, and then review three recent works with high relevance. A picture of the Twistcar toy vehicle appears in Fig. 1a. A child can sit on the vehicle and generate propulsion on level ground by simply applying periodic oscillations of the steering handlebar. A simple planar model of the Twistcar is shown in Fig. 1b, including all notations. The model consists of two rigid Its two-link model with notation links supported by wheeled axles, and connected by a rotary joint. Each link has mass m i and length l i , for i = 1, 2. Each link has moment of inertia J i about its center-of-mass c i (COM), which is located at a distance b i from the wheels' axle. We attach a body-fixed reference frame to each link, such that e i is aligned with the link's longitudinal axis while e ⊥ i is perpendicular to it and parallel to link's wheeled axle. Midpoints on the two wheeled axles are denoted by p i . The relative angle between the links is denoted by φ and the absolute orientation angle of link 1 is denoted by θ . While a mechanical torque τ is applied at the joint through the steering handlebar, it is more convenient to assume that the joint angle φ is directly prescribed and oscillates in time t as φ(t) = φ 0 + ε sin(ωt).

Asymptotic analysis of point-mass model
The first recent relevant work is that of Chakon and Or [29]. They studied a point-mass model of the Twistcar where m 2 = 0 and J 1 = J 2 = 0, and also assumed zero skidding (lateral slippage) of the wheels, which implies where upper dot represents time-derivative d/dt. We now review main observations from [29]. Figure 2 shows typical simulation results from [29] under parameter values of l 1 = 0.5m, l 2 = 0.1m, b 1 = 0.3m, φ 0 = π , ε = 0.6rad, and ω = 1rad/s. Time plot of the speed of rear axle midpoint v 1 (t) =ṗ 1 · e 1 appears as solid line in Fig. 2a. It can be seen that the mean speed has negative mean acceleration (i.e., moving backward), with superposed oscillations of slowly diverging amplitude. Time plot of the vehicle's orientation angle θ(t) appears in Fig. 2b, showing that θ(t) also undergoes oscillations of diverging amplitude. This implies that for long-time simulations, the vehicle's trajectories in (x, y) plane begin to curve into flowerlike loops, as also observed in [40]. Another important effect that has been observed and analyzed in [29] is reversal in direction of the vehicle's net motion depending on the vehicle's structural parameters. Specifically, under oscillations about the straightened configuration φ 0 = π , the direction of net motion depends on the ratio of the steering link's length l 2 to the distance b 1 of the body's COM from the rear axle (see Fig. 1b). For b 1 < l 2 , the vehicle moves in backward direction, whereas for b 1 > l 2 it moves in forward direction, see illustration in Fig. 2c. As an example, the dashed curve overlaid in Fig. 2a shows time plot of v 1 (t) for the same parameter values except for b 1 = 0.05m < l 2 , which gives motion in reverse direction (positive). This nonintuitive phenomenon has been analyzed in [29] using leading-order asymptotic expansion of the point-mass model.

Motion measurement experiments
The second recent work (unpublished) is a series of motion experiments we have conducted with a robotic prototype of the Twistcar vehicle, see Fig. 3a . We used modular parts of Actobotics series from ServoCity in order to construct the robot's structure. The joint was actuated by servo motor powered by on-board lithium-ion battery, and controlled by Arduino Mega controller. We prescribed periodic input command of steering angle φ(t) as in (1), with φ 0 = 0, ε = π/6 and ω = π rad/s. Reflective markers have been attached to the vehicle's links, and the motion has been tracked by an array of infrared cameras using VICON system. Figure 3b shows measurements of the vehicle's forward speed v 1 (t) as a function of time, showing that it first accelerates from rest, and then seems to converge to a periodic solution with bounded oscillations about a mean value. In addition, Fig. 3c  attributed to some small calibration error in the steering angle's servo motor and encoder, as well as asymmetry in the wheels' rotational friction. Nevertheless, the results stand in clear disagreement with the theoretical prediction in [29] of indefinitely increasing mean value of v 1 (t), in addition to diverging oscillations of both v 1 (t) and θ(t), see Fig. 2a and b.
Next, Fig. 3d shows time plots of the measured skid velocities of the rear and front wheel's axles, σ 1 (t) and σ 2 (t). Remarkably, one can see non-negligible skidding, with velocity magnitudes ranging up to 20% of forward mean speed. This clearly indicates that the assumption of ideal no-skid constraints (2) is completely unrealistic and should be relaxed. Moreover, while σ 2 (t) shows alternating reversals in skidding direction, σ 1 (t) seems to display a more complex behavior of transitions between skid and no-skid states. This suggests that a Coulomb-like model of dry friction should be considered. Finally, in experiments where φ(t) oscillated about straightened configuration φ 0 = π , we were unable to reproduce the phenomenon of direction reversal which was predicted in [29] by changing the vehicle's structure and mass distribution. This suggests that the point-mass model is too simplistic, and one should account for the mass and momentof-inertia of both links.

Analysis of the Landshark model
The third recent relevant work is of Bazzi et al [33]. They considered a similar two-link model dubbed as "Landshark," where the two links have non-negligible masses and moments of inertia, such that J i = m i l 2 i and b i = 0. Varying the mass ratio m 2 /m 1 and length ratio l 2 /l 1 , they have identified conditions guaranteeing that the vehicle's momentum grows in a sign-definite direction for any input of steering angle φ(t). In cases where these conditions are not satisfied, they have demonstrated a numerical simulation example where the direction of motion can be reversed for a vehi-cle with specific values of m i , l i by only modulating the input's amplitude and frequency of the steering angle φ(t) oscillating about φ 0 = π as in (1). This interesting effect has not been captured by the leadingorder asymptotic analysis in [29]. Finally, the work [33] also considered a relaxation of the no-skid assumption by using the "skid-angle" model [32,34,35] , which accounts for non-dissipative smooth changes in the skidding velocities. However, this model still predicts unbounded divergence of the forward speed and orientation angle, which again disagrees with our experimental observations in Fig. 3b.

Research goal and statement of contributions
The goal of this work is to extend the theoretical model of the Twistcar vehicle and its analysis in order to explain the discrepancies and disagreement between previous theoretical works and the experimental measurements. More specifically, we make four main contributions which are listed as follows. First, we extend the leading-order asymptotic analysis of [29] in order to incorporate the effect of links' inertias on direction of net motion and its reversal. Second, we derive the next-order expansion term of the vehicles' speed in order to find necessary and sufficient conditions for achieving direction reversal by changing only the oscillation amplitude ε of the steering angle. Third, we generalize the theoretical model in order to incorporate Coulomb friction bounds and hybrid dynamics of skidstate transitions, and show numerical simulation results that qualitatively agree with our experimental measurements. Finally, we numerically analyze the influence of input frequency on characteristics of periodic motion with skid-state transitions, and obtain optimal frequencies that maximize mean speed as well as energetic efficiency.
The rest of the work is organized as follows. The next section formulates of the Twistcar's dynamic equations as a nonholonomic system under ideal no-skid constraints, and presents some numerical simulation results. Section 3 conducts asymptotic analysis including links' inertias and next-order expansion, and highlights the effect of links' inertias and input's amplitude on direction reversal. Section 4 provides expressions for normal reaction forces assuming "flat plane" model of the vehicle, and introduces Coulomb's friction constraints and the hybrid dynamics of skid-state transi-tions. Section 5 presents numerical analysis of periodic motion and the influence of input's frequency. Finally, Section 6 concludes the work.

Problem formulation with ideal constraints
We now present formulation of the system's nonholonomic dynamics under ideal no-skid constraints. We choose generalized coordinates q = (x, y, θ, φ) T , where x, y denote the position c 1 of the COM of link 1, θ is the link's orientation angle, and φ is the steering joint angle, see Fig. 1b. The skidding velocities σ 1 , σ 2 defined in (2) can be expressed as Defining the matrix W(q) = [w 1 (q) w 2 (q)] T , the nonholonomic no-skid constraints can be written as The system's constrained dynamic equations of motion can be obtained using Euler-Lagrange derivation as (cf. [41]): where L = T − U is the system's Lagrangian, T is the total kinetic energy, U is the potential energy, F q is the vector of generalized forces and = (λ 1 , λ 2 ) T is the vector of constraint forces (Lagrange multipliers).
Since the vehicle here moves in horizontal plane, the system is not governed by any changes in potential energy and U can thus be ignored. The kinetic energy for planar motion of the 2-link system can be obtained as where full expressions for elements of the matrix of inertia M(q) are given below in (8). Using (5) and (6), the equations of motion can be written in a standard matrix form as and τ (t) is the actuation torque at the steering joint. Explicit expressions for the elements of M(q) and the vector of velocity-dependent terms B(q,q) from (7) are given as: Differentiating the velocity constraints in (4) with respect to time gives: Combining equations (7) and (9) gives a differential algebraic system where the vector of constraint forces can be eliminated, and thus the accelerationsq can be integrated numerically under given input torque τ (t). Next, we consider the case where the steering angle is assumed to be a directly controlled input, so that φ(t) and all its time-derivatives are prescribed functions. We note that the torque τ appears only in the fourth row of (7). Thus, this row can be ignored, so that (7) and (9) now give a 5 × 5 linear system in the unknowns {ẍ,ÿ,θ, λ 1 , λ 2 }, which can be numerically integrated to solve for x(t), y(t) and θ(t). Note that one may choose to formulate the system dynamics using Appell-Gibbs approach [42,43], which gives a more elegant elimination of the constraints [13,14]. Nevertheless, we choose here to use Lagrangian formulation, following [29], since it enables explicit formulation of constraint forces which are crucial for checking friction inequalities, as detailed later We now present numerical simulation results of the Twistcar. Physical parameter values which correspond to the robotic vehicle in Fig. 3a are chosen as  It can be seen that the condition b 1 < l 2 derived in [29] is no longer sufficient for reversing the direction of net motion. However, changing the mass and inertia of link 1 to m 1 = 4.3, J 1 = 0.025K g · m 2 (equivalent to adding a uniform disk of mass 2K g and radius 0.03m placed at the link's COM), leads to direction reversal as shown in the dash-dotted curve of v 1 (t) in Fig. 4a. Analytical derivation of exact conditions for direction reversal with masses and inertias is presented in the next section. From Fig. 4b, it can be seen that the constraint forces λ i (t) are diverging in time. Moreover, since λ i (t) scale as ω 2 with actuation frequency, fast oscillations of the input φ(t) require large constraint forces which may exceed friction limitations, practically resulting in wheels skidding. This issue is addressed in Sect. 4. Finally, oscillations of the steering torque τ (t) are also diverging in time, which in unrealistic in any actuated robot due to saturation bounds on the motor's torque.

Reduction and asymptotic analysis
In this section, we conduct reduction of the Twistcar's dynamic equations of motion (7) and then present asymptotic analysis under small oscillation amplitude ε of the steering angle input φ(t). This is a generalization of the analysis presented in [29] for the point-mass model.

Reduction in the dynamic equations
First, we define nondimensional parameters: The system's velocities, accelerations, forces and torques are then normalized by characteristic scales of The vector of generalized body velocities is then defined as (7),(9) and pre-multiplying (7) by R T , one obtains the body-frame equations of motion as: and we used R T E = E. Note that this process eliminates the dependence of (12) on the orientation angle θ due to invariance of the system with respect to rigid-body transformation. The expressions for matrices from (12) are given as: The last reduction step is elimination of the constraints and constraint forces. This is done by defining a vector of constraint-free velocities v f = (u v) T as: Note that the vectors w u (φ), w v (φ) are mutually orthogonal and satisfy the constraints, that is, The constraint-free velocities in (13) have the same role as pseudo-velocities in Appell-Gibbs formulation, cf. [14]. The physical meaning of the constraint-free vector field w v (φ) is motion along an arc while the steering angle remains constant (φ = 0). The second vector field w u (φ) is orthogonal to w v and also satisfies the constraints. That is, it involves varying steering angle such that u is proportional toφ up to the scaling factor g(φ) given in (14). Substituting (13) into (12) and pre-multiplying by S T gives the reduced equations: Note that the constraint forces have been eliminated from (15). Assuming that the steering angle φ(t) is a prescribed function of time, u(t) can be extracted from the relation (14). Moreover, the structure of S and E implies that the steering torque τ does not appear at the second row in equation (15). This gives a scalar firstorder differential equations that governs the dynamic evolution of the generalized unconstrained velocity v, taking the form: The function F in (16) is highly complicated and its explicit expression is not shown here. It depends quadratically onφ, linearly onφ, v, and trigonometrically on φ. Moreover, F is even in φ(t), which implies that: Similar symmetry relation also holds about φ = π .

Asymptotic analysis-perturbation expansion
We now derive approximate solution of (16) via perturbation expansion, assuming small amplitude ε of input oscillations. We expand v(t) as a power series: . Next, we expand F(·) in (16) as The dependence on the angle φ is expanded about φ 0 ∈ {0, π}. The expansions ofv and F from (18) and (19) are substituted into (16), while φ(t) and its derivatives are taken from (1) under time-scaling ω = 1. We have used symbolic computations in MATLAB for this process. The resulting equation is then rearranged by different powers ε k in order to obtain a sequence of differential equations for v k (t) which can be solved iteratively, as follows. The zero-order term givesv 0 = 0; hence, we obtain v 0 (t) = v 0 = const. Due to the even-function symmetry of F(·) in (17), all odd powers also givev k = 0 for odd k; hence, these terms are neglected. The ε 2 -terms give the leading-order dynamics of the formv 2 (t) = F 2 (v 0 , t). Integrating in time (assuming zero initial conditions) gives: where a 2 , b 2 are constant coefficients that depend on the system's parameters, the value of φ 0 = {0, π} and v 0 . Explicit expressions for the coefficients are given in Table 1, where we assume starting from rest, v 0 = 0, for simplicity. Note that a 2 reduces to the simple expression in [29] for the point-mass case, by substituting κ = η 1 = η 2 = 0. Expansion to the next order of ε 4 -terms, which has not been considered previously in [29], gives equation of the formv 4 (t) = F 4 (v 0 , v 2 , t).
Assuming v 0 = 0, integrating in time and substituting the expression for v 2 from (20), one obtains: The expression for a 4 appears in Table 1, whereas the lengthy expressions for b 4 , d 4 (associated with higherorder oscillating terms) are not shown. Using (20) and (21), the expression for v(t) has oscillating parts superimposed on a constant-rate accelerating mean partv(t), which is expanded up to order 4 as The mean accelerationā = dv/dt of the vehicle is thus given as: The influence of system's parameters and input amplitude ε onā, and particularly on its sign reversal, is analyzed next.

Small-amplitude parameter-dependent direction reversal
Assuming small amplitude ε 1 of the oscillating steering angle φ(t), the leading-order term of mean acceleration in (23) is a 2 . For oscillating about φ 0 = 0, substituting ζ = 1 into the expression for a 2 in Table 1, it is easy to see that it is always positive. However, oscillating about φ 0 = π (i.e., ζ = −1) gives possible conditions for sign reversal of a 2 depending on the vehicle's parameters. For simplicity, we assume that the COM of link 2 is located at its midpoint, β 2 = 0.5. Substituting this and ζ = −1, we obtain: As an example, Fig. 5a plots the leading-order acceleration a 2 as a function of β for α = 0.5 and κ = 0 in three solid-colored curves corresponding to three different choices of other parameter values. Case (i) is "point-mass" model η 1 = κ = 0 as studied previously in [29]. It can be seen that the direction of net motion reverses to a 2 > 0 for β 1 < α. Case (ii) adds mass and uniform inertia of the links η 1 = η 2 = 1/12, κ = 0.22 as in the experimental robot (Fig. 3a). It can be seen that the effect of added inertias precludes the possibility of direction reversal (a 2 < 0 for all β 1 ). This explain the difficulty to achieve direction reversal in  our experiments. Case (iii) corresponds to adding a uniform disk to link 1 as described previously in the numerical simulations (Fig. 4a) which corresponds to η 1 = 1/22, κ = 0.1. In this case, there is a narrow range of β 1 for which the sign of a 2 is reversed. All this can be explained by the leading-order expression in (24), which highlights the separate contribution of each parameter of links' inertia. Finally, the dashed colored curves denote mean accelerationā/ε 2 which is computed via numerical simulations for input amplitude of ε = 0.6rad, showing good agreement with the leading-order expression in (24). (Mean acceleration was numerically computed by taking the integrated v(t) and smoothing it using a moving-average filter, and then applying linear regression to the resulting smoothed signal in order to calculate mean slope).

Amplitude-dependent direction reversal
In case where the input's amplitude ε is not very small, the vehicle's mean accelerationā depends on both a 2 and a 4 , as given by (23). When they differ in signs, i.e., satisfy a 2 a 4 < 0, for small amplitudes ε, the sign ofā is determined by that of a 2 , whereas for increased amplitudes, the influence of a 4 causes sign reversal. The critical amplitude ε 0 for sign reversal can be obtained from (23) as Moreover, for small amplitudes ε < ε 0 , it is easy to show that there exists an optimal amplitude ε * = ε 0 / √ 2 for which |ā| attains maximal value. We now revisit the "Landshark" model and analyze the amplitude-dependent direction reversal observed in [33]. Using the definitions in (10), the Landshark's nondimensional parameters are η 1 = η 2 = 1, β 1 = β 2 = 0, leaving only the links' length and mass ratio α, κ as variables. Substituting this and ζ = −1 (i.e., oscillations about φ 0 = π ) into the expressions for a 2 and a 4 in Table 1 gives: Figure 5b plots regions of ±a 2 , ±a 4 > 0 in the plane of κ, α. The shaded area denotes region of a 2 a 4 < 0, where amplitude-dependent direction reversal can occur. The result extends the analysis in [33], which only provided necessary conditions for possible direction reversal. For comparison, the hatched regions in Fig. 5b were labelled as '±' in [33], and satisfy necessary conditions for sign reversal under some input of φ(t). Figure 5b shows that these regions include the shaded regions satisfying our necessary and sufficient condition a 2 a 4 < 0.
As specific examples, we fix κ = 1.8 and choose three values α = {0.5, 0.57, 0.75}, which are marked by '×, +, * ' in Fig. 5b, respectively. For each value of α, we conduct numerical simulations with varying amplitudes ε and compare the numerically computed mean accelerationā (dashed curves) to its four-order approximation in (23) (solid curves), all plotted as a function of ε in Fig. 5c. It can be seen that only the value corresponding to marker "+" which lies within the shaded region in Fig. 5b exhibits direction reversal as a function of amplitude ε in Fig. 5c. In contrast, the values corresponding to markers "×, * " which lie {inside, outside} the hatched region of [33] in Fig. 5b do not result in direction reversal. Finally, Fig. 5c shows that the deviations between numerical and 4 th -order asymptotic approximation (dashed and solid curves) are vanishing for small ε and begin to increase for ε > 1, as expected. This concludes the first part of this work, which assumes ideal no-slip constraints. In the next section, we relax this assumption using Coulomb's friction.

Friction-bounded hybrid dynamics of skid-state transitions
We now present bounds on the constraint forces of the wheel-ground contact according to Coulomb's friction model, and then formulate of the vehicle's hybrid dynamics under friction-induced skid state transitions.

Coulomb's friction constraints and normal forces
For each wheel's axle, the constraint force λ i has to satisfy Coulomb's friction bound given by |λ i | ≤ μN i (27) where N i is the normal reaction force acting at the axle's wheel(s) inẑ direction, and μ is Coulomb's coefficient of dry friction. Positions of the left and right wheels on the rear axle are denoted by p 1a , p 1b , see Fig. 1b. The two normal forces acting on those wheels are denoted as N 1a , N 1b , such that the total normal force on the back axle is For the front axle, the normal force N 2 acts on the single wheel positioned at p 2 . We now make a simplifying assumption that the vehicle is flat and the links lie within the horizontal x y plane with z = 0. This assumption holds when the actual height of the links' centers-of-mass is much smaller than the lengths l i . Under this assumption, accelerations of the links in horizontal plane generate negligible moments about the wheels-ground contact points. This implies that the three normal forces {N 1a , N 1b , N 2 } are decoupled from the vehicle's dynamics in x y plane and can be obtained from statically determinate equilibrium equations of total {force, torque} balance {along z, about x y } directions, which are written as: This is a linear 3 × 3 system in the unknowns {N 1a , N 1b , N 2 }, which is expressed in body-fixed frame as: The solution of (29) for normal forces depends only on the steering angle φ. A singular case of unbounded forces occurs when l 1 = l 2 cos φ, but this case is ruled out assuming l 1 > l 2 . In addition, tipover of the vehicle due to vanishing of one normal force may also occur, characterized by the two links' common COM lying on the line connecting two wheels. This case is also rare for reasonable values of the vehicle's geometric structure (more details appear in the thesis [44]).

Hybrid dynamics with skid-state transitions
We now formulate the system's hybrid dynamics which accounts for friction bounds and transitions between skid state. First, consider the case where the no-skid constraint (3) at the i th wheel axle is satisfied. The constraint force λ i must satisfy Coulomb's dry friction inequality (27), where the normal reaction force N i (φ) is obtained from solving (29). When the friction bound in (27) is reached for the i th axle, it begins to skid, so that σ i (t) = w i (q) ·q = 0, and the constraint force reaches its maximal value, λ i = −μ sgn(σ i )N i . The skid states described above lead to 2 2 = 4 possible "modes" of the system-combinations of the axles' skid states. The dynamics of each mode is governed by a different set of equations, described as follows. First, Lagrange's equations (7) still hold, giving three scalar equations (assuming φ(t) is prescribed and eliminating the third row in (7) which also contains the torque τ (t)). However, the equation (9) of the time-derivative of the no-skid constraints does not always hold, and should be replaced by two equations which depend on the skid state of each axle, given as: axle i is skidding (30) for i = 1, 2. Combining equations (30) with (7) gives a system of differential-algebraic equations which yields a 5 × 5 linear system in the unknowns, {ẍ,ÿ,θ, λ 1 , λ 2 } for a prescribed angle input φ(t).
Transitions between different modes of the hybrid dynamics are described as follows. A non-skidding axle begins to skid whenever its constraint forces reaches its frictional bound, |λ i | = μN i . When the slip velocity of a skidding axle vanishes w i (q) ·q = 0, the wheel may switch back to no-skid state, provided that the constraint force at the initial instant right after switching satisfies its bound |λ i | ≤ μN i . Otherwise, the wheel's slip velocity reverses its sign. All this gives rise to a hybrid dynamical system, containing transitions between different skid states.
Singularity of the hybrid dynamic equations should also be examined, by constructing the matrix of the linear system which governs the dynamics under each mode of skid states, and checking whether its determinant can cross zero. The details of this singularity analysis appear in Appendix, for both cases of controlled actuation torque τ (t) or controlled steering angle φ(t). Summarizing this analysis, it is concluded that for each state that involves skidding, the simplifying assumptions of point-mass model, κ = η 1 = 0 from [29], may lead to singularity. This justifies considering the two links' masses and inertia in our work.

Numerical simulation results of the hybrid dynamics
We now present numerical simulation results of the hybrid system. The simulations are conducted using ode45 procedure of MATLAB for adaptive-step integration, with event functions for detecting zerocrossing conditions for each mode's termination. Parameter values are chosen based on the robot prototype in Fig. 3a, as m 1 = 2.29K g, m 2 = 0.229K g, l 1 = 0.362m, l 2 = 0.19m, b 1 = 0.1249, b 2 = 0.5l 2 , I i = m i l 2 i /12 for i = 1, 2, d = 0.1m and friction coefficient of μ = 0.3.

Periodic solutions with skid-state transitions
First, we show how the system's motion converges to periodic solutions with skid-state transitions. The input steering angle is φ(t) = ε cos(ωt) with ε = 0.6rad and ω = 2.5rad/s, under initial conditions of zero velocityq(0) = 0. Figure 6a,b shows time plots of the forward speed v 1 (t) and orientation angle θ(t). Pieces of dashed curves represent motion, while the front axle is skidding, σ 2 (t) = 0 whereas solid curves correspond to no-skid motion. It can be seen that the system begins with no-skid motion and accelerates until onset of skidding. Then, the solution converges to periodic motion with skid-state transitions, which contains bounded oscillations of v 1 (t) and θ(t) about mean speed and orientation, respectively. This result is remarkably different from the unrealistic solution divergence of the no-skid model [29], see Fig. 2a,b. Figure 6c shows time plots of the forces' ratio λ i /(μN i ) vs. time for i = 1, 2, overlaid with the skid velocities σ i (t). It can be seen that the constraint forces converge to periodic solution while staying bounded by the friction inequalities (27), in contrast to their divergence in the no-skid simulation (Fig. 4b).
The slip velocity of the front wheel σ 2 (t) is initially zero, and then becomes nonzero whenever the constraint forces reaches its bound |λ 2 (t)| = μN 2 (t). The motion converges to periodic solution, while the rear axle never skids, σ 1 (t) = 0, since the associated constraint force is always below its friction bound, |λ 1 (t)| < μN 1 (t). One can see that in steady state, the angle θ(t), contact forces λ i (t) and slip velocity σ 2 (t) are periodic in t p = 2π/ω. Moreover, they satisfy half-period anti-symmetry Symmetry of the vehicle about φ = 0 implies that the forward speed is doubly periodic v 1 (t + t p /2) = v 1 (t). This also agrees with the asymptotic analysis of the no-slip system in [29]. (See also (20) where the frequency of forward speed is doubled).
Next, we consider simulation results with higher actuation frequency of ω = 6rad/s. The system's motion converges to a steady-state periodic solution with significantly different characteristics, as shown in the time plots of Fig. 7. One can see that the rear axle now undergoes alternating sequence of stickslip transitions where σ 1 (t) switches between zero and nonzero values . On the contrary, the front axle always slips, σ 2 (t) = 0, while its direction reverses alternately, every half period. This type of motion is remarkably similar (qualitatively) to the experimental measurements obtained with the robotic vehicle, see Fig. 3. Upon changing the actuation frequency, the mean value of v 1 (t) displays a small change (Figs. 6a, 7a). On the other hand the oscillation amplitude of the vehicle's orientation angle θ(t) changes drastically Figs. 6b, 7b, respectively, where θ = θ max − θ min ). In both cases, the angle θ(t) oscillates about a constant mean value, which implies net translation along a straight line with zero net rotation .

Influence of input frequency on performance
We now analyze the influence of input frequency ω on properties of the hybrid periodic solution. We conduct numerical simulations under input (1) with φ 0 = 0, ε = 0.6rad and frequency ω that varies in small incre-  Figure 8a plots the percentage that each skid state occupies out of the period time t p = 2π/ω. It can be seen that for low frequencies, the only existing states are no-skid (σ 1 = σ 2 = 0, solid curve) and stick-skid (σ 1 = 0, σ 2 = 0, dashed curve). As the frequency ω increases, the no-skid state become shorter. At ω ≈ 3.25rad/s, the state skid-skid emerges (σ 1 , σ 2 = 0, dotted curve), and its relative time duration grows with ω. At ω ≈ 3.5rad/s the no-skid state disappears. The "noisy" values around ω = 3.5 are a result of numerical sensitivity of the event-based integration when two different events are almost coinciding. The vertical dash-dotted lines in Figs. 8, 9 represent the frequencies ω = {2.5, 6} from the simulations in Figs. 6, 7, respectively.
In the steady-state periodic solution, the body orientation angle θ(t) oscillates about a mean valueθ (see Figs. 6b, 7b), so that the net motion is pure translation with zero net rotation. Figure 8b plots the oscillation amplitude of θ(t) as a function of frequency ω, showing that it decays with ω. Next, we define the forward displacement per cycle as S = c 1 (t + t p ) − c 1 (t) ·ē 1 , whereē 1 = (cosθ, sinθ) T is unit vector along the mean body orientation. The mean translation speed is defined as v = S/t p . Plots of the displacement S (solid) and speed v (dashed) as a function of frequency ω are overlaid in Fig. 9a. One can see direction reversals for low frequencies, caused by large oscillation amplitudes of θ(t) which result in large rotations θ of more than a full revolution of the body orientation during cycle for ω < 0.9rad/s. These large rotations are not physically realistic and did not occur in the experiments (Fig. 2), as they require very small amount of energy dissipation due to skid during the cycle. In addition, one can see in Fig. 9a extremum points of optimal frequencies for which |S| and |v| are maximized. Finally, we compute the mechanical energy expenditure per cycle W = τ (t)φ(t)dt. A plot the energetic cost of transport W/|S| as a function of frequency appears in Fig. 9b. One can see three minimum points of optimal energy efficiency, corresponding to the three regions of reversed directions. If one considers only the more practical range of sufficiently large frequency with reasonable rotations, there are unique optimal frequencies for maximizing distance, speed and energy efficiency. (More precisely, the frequency ω should be replaced by the nondimensonal quantity ψ = lω 2 μg which represents ratio of characteristic constraint forces to normal reaction forces. The ratio ψ actually governs the system's dynamic behavior. That is, increasing ω 2 is equivalent to decreasing the friction coefficient μ). An open question which requires further exploration is whether this system has co-existing hybrid periodic solutions, stability transitions, bifurcations and chaos, cf. [45]. In our analysis, we focused only on periodic solutions that also possess half-period symmetry and involve zero net rotations. In such solutions, we did not find any co-existing solutions or stability loss of our solutions, which was assesed by numerical evaluation of Poincaré map linearization eigenvlaues. It could be possible that periodic solutions with nonzero net rotations exist, and the search for them and their stability analysis and bifurcations are left for future investigation.

Conclusion
We have presented significant extensions of the theoretical model of the Twistcar's nonholonomic dynamics under oscillating steering angle. The point-mass model from our previous work [29] has been extended to include mass and inertia of both links, and the asymptotic expansion enabled studying the effect of the vehicle's geometry and inertia parameters on reversal in direction of net motion. In addition, carrying the asymptotic expansion to next-order correction terms enabled discovering cases of direction reversal depend-ing on the input amplitude ε, which provides explanation to results from [33] for the Landshark model. Next, we have relaxed the no-skid assumption of ideal nonholonomic constraints by considering Coulomb friction bounds on constraint forces, and introducing hybrid dynamics of skid-state transitions. Numerical simulations of the complex hybrid dynamics show convergence to periodic solution with skid-state transitions and bounded oscillations of both forward speed and body orientation angle, in agreement with previous experimental observations. Finally, the influence of actuation frequency ω on the vehicle's performance in steady-state periodic solution was studied. We have found very close optimal values of frequencies that maximize the mean speed, travel distance per cycle, and energetic cost of transport (energy per distance).
We now briefly discuss limitations of the current results and sketch possible directions for future extension of the research. First, while the theoretical predictions agree qualitatively with observations from our preliminary experiments, a more thorough quantitative analysis of controlled motion experiments is crucially needed in order to study effects of actuation frequency and amplitude as well as vehicle's structural parameters on its motion. Second, general models of frictional contact forces and other possible mechanisms of energy dissipation such as rolling resistance should be studied both experimentally and in mathematical analysis. From our preliminary experiments, the value of Coulomb's dry friction coefficient is crudely estimated within the range 0.1 ≤ μ ≤ 0.7. However, it depends crucially on material properties of the terrain and the wheels' structure, and may contain directional anisotropy. Thus, a more thorough analysis of the contact mechanics and dissipation is needed. In addition, the simplifying assumption of "flat vehicle" with negligible height of its COM should be relaxed in order to investigate the influence of COM accelerations on ground contact forces in dynamic maneuvers. Third, extending the small-amplitude asymptotic analysis from Section 3 and from [29] to the case of hybrid skid-states motion is a highly challenging analytical problem. Fourth, the current analysis is limited to the case of harmonic input of steering angle. Finding optimal gaits of arbitrary time-periodic inputs is still an open problem for this type of system. One may consider a purely computational approach as in [46], where the input is taken as a truncated Fourier series whose finite set of coefficients is optimized for maximal dis-placement. Alternatively, one can apply optimal control formulation as in [25] or its direct numerical solution as in [47]. Finally, the research should be extended to account for multi-link wheeled vehicles, possibly with passive elastic joints, as in [48][49][50]. Availability of data and material Not applicable.

Conflict of interest
The authors have no conflict of interest.