Energy-conserving intermittent-contact motion in complex models

The possibility of a persistent motion, such as hopping or walking on a level ground, in a mechanical system with inelastic collisions may appear paradoxical, at first glance. However, no laws of mechanics are violated if all the intermittent contacts along the trajectory occur at a vanishing (at the point of contact) velocity. For a solution with purely repulsive ground reaction force, as in conventional walking or hopping, acceleration must also vanish at the contact [1]. We call such motion collisionless [2].

Although collisionless solutions have been considered in one- [3] two- [2, 4] and three- [5] dimensional models, analytical consideration was effectively limited to 2-DOF complexity. For example, the three-dimensional bipedal walker in [5], while nominally being a 5-DOF model, consisted only of two articulated parts, conveniently decoupling into 2-DOF subsystems in the harmonic approximation (plus a 1-DOF subsystem with a forcing term). In the analysis of those simple models it was noted that although a complete description of a harmonic system includes both normal frequencies and normal vectors, the core equations of the collisionless motion problem could be formulated exclusively in terms of the normal frequencies [3, 5]. Whether this peculiarity is a general property, not limited to low complexity models, remained unknown.

Not every model admits a collisionless solution for any set of model parameters. However, a counting argument suggests that the dimensionality of constraints (needed to meet collisionless solution requirements) on the model parameter values is small and is independent of the model complexity [5].

On the one hand, this may indicate the relevance of collisionless models for understanding animal locomotion [6], as well as for designing efficient walking machines [7]. Indeed, the mechanical cost of transport (COT) – the standard measure of locomotion efficiency equal to energy expended per weight per distance traveled – is remarkably low for humans, dogs and other animals at walking speeds (0.08 for humans and 0.04 for dogs, according to [8]). In contrast, the majority of man-made bipeds and quadrupeds today (with some notable special purpose built exceptions [9, 10]) typically demonstrate an order of magnitude higher COT.

On the other hand, looking numerically for a suitable collisionless solution, in even a low-DOF non-linear model, is not a simple task [4]. And the task may be much harder for higher-DOF models. Therefore, it is desirable to extend the analytical treatment to higher complexity models, potentially gaining useful insights along the way.

In this work, we extend the analytical analysis to $N$-DOF linear (i.e. in harmonic approximation) models. We show that the formulation in terms of the normal frequencies is a general property of collisionless motion, not limited to low complexity models. Starting with $N=2$ example, we conjecture a general collisionless solution existence condition. Compelling evidence in support of its correctness is then provided by analyzing a general $N$ solution in the critical region, where the collisionless solution ceases to exist. We illustrate how the search for the roots of collisionless equations can be aided by two-dimensional plots. The practical value of the harmonic approximation solution is that it can be straightforwardly evolved numerically toward a collisionless solution of a non-linear model [5].

The paper is organized as follows. In Section 2 we more rigorously outline the concept of collisionless motion, problem specifications and related notations. The main equations are derived in Section 3 for general $N$. Known solutions for $N=2$ are reproduced in Section 4. In Section 5 the critical region solution is analyzed. The developed approach is applied to solve a new $N=3$ model in Section 6. We conclude with an outlook for further research in Section 7.

We consider a motion of a mechanical system interacting with the ground via purely inelastic collisions and non-slipping contacts. We assume no internal collisions in the system, collisions may only occur between its parts and the ground. When a collision happens at a finite velocity, a finite amount of mechanical energy dissipates at the collision. In the rest of the paper the term collision will only be used in this context. If, on the other hand, velocity continuously vanishes at the moment of contact, so no energy gets lost, we call it an impact. Unlike a collision, an impact is time-reversible. We refer to the motion containing impacts and no collisions as collisionless.

A collisionless motion can be viewed as a temporal sequence of phases, each phase characterized by a distinct configuration of contacts and separated from the adjacent phases by impacts. Contacts constrain the system, reducing its number of DOF. The impact dimensionality is defined as the number of DOF that get constrained at the impact. In this paper we limit the consideration to a periodic collisionless motion with two phases separated by a one-dimensional impact. The phase with lower DOF is called constrained, while the other is called unconstrained. We will use a prime to distinguish quantities in the constrained phase, (e.g. $x$ in the unconstrained phase becomes $x^{\prime}$ in the constrained phase). We will use $p$ superscript to represent a quantity in either of two phases, (i.e. $x^{p}$ can be $x$ or $x^{\prime}$). If a statement is obviously applicable to both phases, we may optionally omit $p$.

Following the literature [3, 2, 5], we restrict the collisionless solution to a particular class of periodic trajectories characterized by two symmetry points $P^{p}$ associated with the constrained and unconstrained phases. At $P^{p}$, the trajectory must be continuous and invariant under a simultaneous time reversal $T$ and a spatial transformation $S^{p}$ which preserves the Euclidean metric and leaves the ground invariant. By the invariance of the Lagrangian to $S^{p}$ and the time-reversibility of collisionless motion, the invariance to $TS^{p}$ extends from $P^{p}$ to the whole trajectory. The imposed symmetries simplify the consideration, as the number of independent variables is reduced and one only needs to consider the trajectory between $P$ and $P^{\prime}$, which can then be unfolded into the full solution according to the symmetries.

Let $x$ be an $N$-dimensional coordinate vector of the unconstrained phase. Let $x_{N}$ be the component of $x$ that gets constrained at the impact, so $x_{N}$ is constant in the constrained phase. To realize a collisionless motion, in addition to $\dot{x}_{N}\to 0$ at the impact, it is also required that $\ddot{x}_{N}\to 0$, otherwise the contact cannot persist [1]. Let $t_{imp}^{p}$ be the time of the impact (laying between the adjacent $P$ and $P^{\prime}$) measured from the symmetry point $P^{p}$, and $\tau^{p}\equiv|t_{imp}^{p}|$. In summary, the following conditions must be satisfied at the impact:

In the rest of the section, we briefly review some matrix-related notations and conventions employed in the paper. See A for more details.

We use columnar format for vectors. In a matrix $A$ and vector $a$: $A_{i}$ is the $i$th matrix row, $a_{i}$ is the $i$th vector element, $A^{i}$ is the $i$th matrix column and $A_{i}^{j}=A_{ij}$ is the matrix element at the intersection of $A_{i}$ and $A^{j}$. We use square brackets to denote a matrix, commas to separate columns and semicolons to separate rows. For example, for a $n\times m$ matrix $A$, $A=[A^{1},A^{2},...,A^{m}]=[A_{1};A_{2};...;A_{n}]$ and $A^{i}=[A_{1i};A_{2i};...;A_{ni}]$. A matrix superscript should not be confused with a power notation, which is only used for scalar quantities, i.e. $A_{ij}^{k}=(A_{ij})^{k}$. We use $A_{(i)}$ and $A^{(j)}$ to denote matrices obtained from $A$ by removal of $A_{i}$ and $A^{j}$ respectively. Also, $A_{(ij)}\equiv A_{(i)}^{(j)}$. For a vector or matrix $a=[a_{1};a_{2};...;a_{N}]$, we define $\bar{a}=a_{(N)}$. We reserve $I$ for $N\times N$ identity matrix and $\mathbbm{1}$ for $N$-dimensional vector of ones.

For matrices $a$ and $b$ we denote the operations of elementwise multiplication as $a\cdot b$ (also known as Hadamard product) and elementwise division as $a/b$. The operations are defined to support broadcasting and have a higher operator precedence than ordinary matrix multiplication, e.g. $ab\cdot c=a(b\cdot c)$. We introduce these notations to reveal a peculiar and somewhat simple structure of the derived equations, that may be not obvious otherwise.

Consider a system of $N$ linear oscillators $x=[x_{1};x_{2};...;x_{N}]$ with the kinetic energy $T$ and potential energy $V$ given by

where $m$ is positive definite and $k$ is non-singular. The corresponding Lagrangian in the presence of an external force $F$ is

The system is invariant under a simultaneous static shift in $x\to x+x^{0}$ and $F\to F+F^{0}$ with $F^{0}=kx^{0}$. One can switch to the basis of normal coordinates $Q$ via a coordinate transformation $x=XQ$ to find

where the columns of $X$ are the eigenvectors of $m^{-1}k$ with the corresponding eigenvalues forming the spectrum vector $\lambda$, the force $f=X^{\mathsf{T}}F$, and the constant $c$ is selected such that $X_{N}^{N}=1$. We assume that all $\lambda_{i}$ are different (i.e. the spectrum is non-degenerate) and arranged in ascending order, for convenience. Each $\lambda_{i}$ corresponds to a normal mode oscillating with a frequency $\omega_{i}=\sqrt{\lambda_{i}}$. When $\lambda_{i}<0$, the motion is unbounded, diverging with time as $\exp(\pm\nu_{i}t)$, where $\nu_{i}\equiv|\omega_{i}|$.

As can be easily verified from the equations of motion following from the Lagrangian in Eq.(4), given a force $F(t)=F^{\upsilon}(t)=F^{\upsilon}(0)e^{i\upsilon t}$ oscillating with a frequency $\upsilon$, the forced oscillations $Q^{\upsilon}$ are

where $f^{\upsilon}=X^{\mathsf{T}}F^{\upsilon}$. In the original coordinates, $x^{\upsilon}=XQ^{\upsilon}$. Consider now a force applied only to $x_{N}$, oscillating with a normal mode frequency $\omega^{\prime}_{i}$ of the constrained system, (that is $F_{i\neq N}(t)=0$, $F_{N}(t)=F^{\omega^{\prime}_{i}}_{N}e^{i\omega^{\prime}_{i}t}$ and $\omega^{\prime 2}_{i}=\lambda^{\prime}_{i}$). Then $x^{\omega^{\prime}_{i}}_{N}=0$, while $\bar{x}^{\omega^{\prime}_{i}}$ is the normal mode of the constrained system corresponding to $\lambda^{\prime}_{i}$. Note that from Cauchy’s interlace theorem [11] it follows:

Let $\tilde{x}^{\omega^{\prime}}=x^{\omega^{\prime}}/cF_{N}^{\omega^{\prime}}$ and $X^{\prime}=[\tilde{x}^{\omega^{\prime}_{1}},\tilde{x}^{\omega^{\prime}_{2}},...,\tilde{x}^{\omega^{\prime}_{N-1}}]=X\cdot X_{N}(1/(\lambda\bar{\mathbbm{1}}^{\mathsf{T}}-\mathbbm{1}\lambda^{\prime\mathsf{T}}))$. Then the transformation between the normal coordinates of the constrained system¹¹1$X^{\prime}$ does not follow the same normalization convention as $X$, so the Lagrangian of the constrained system does not have the form of Eq.(4), when written in $Q^{\prime}$. and the original coordinates is $x=X^{\prime}Q^{\prime}$, while $X^{\prime}_{N}=0$. If in the above analysis the constraining force also has a static component $F_{N}^{0}$ (that is $F_{N}(t)\to F_{N}(t)+F^{0}_{N}$), then the transformation becomes $x=X^{\prime}Q^{\prime}+x^{0}$, where $x^{0}=(k^{-1})^{N}F^{0}_{N}$. We can write the expression for $X^{\prime}$ and the relation $X^{\prime}_{N}=0$ in terms of a matrix $M$

as

and

Let $\eta^{\mathsf{T}}=X_{N}\cdot X_{N}$. Since $\eta_{N}=X_{NN}^{2}=1$, $\eta$ can be determined from Eq.(9):

Note that $M$ is a so-called Cauchy matrix, and $\bar{M}$ is a square Cauchy matrix, for which determinant and inverse are given by explicit formulas [12], that are evaluated in $\mathcal{O}(N^{2})$ operations.

Normal modes are orthogonal in the sense that their evolution (with time) is governed by the dynamics of a one-dimensional harmonic oscillator. That is, the time dependence of a normal coordinate $Q_{i}$ can be parameterized by two constants $q_{i}$ and $s_{i}$ as $Q_{i}(t)=q_{i}g(t,\omega_{i},s_{i})$, where

$s^{2}\leq 1$ and $s_{\lambda}=\textrm{sign}(\omega^{2})$. A complete solution of the unconstrained and constrained systems, $x(t)=XQ(t)$ and $x^{\prime}(t)=X^{\prime}Q^{\prime}(t)+x^{0}$ respectively, are then readily available, given $q^{p}$ and $s^{p}$ constants. In the context of our search for a collisionless solution, we will see that $s^{p}$ are fixed from the outset by the symmetry of the solution, while $q^{p}$ are determined from equations on $x^{p}$ and their time derivatives at the impact. Interestingly, the equations for the impact times, as well as $q^{p}$, can be expressed through the spectra $\lambda^{p}$ alone, without the use of $X^{p}$. This may look surprising because, for a given $q^{p}$, the configuration $x^{p}(t)$ generally depends on $X^{p}$. In our case, however, this dependence is confined to $\bar{x}^{p}(t)$, as follows from Eq.(10). In other words, $X^{p}$ only influences the part of the configuration that is irrelevant to the impact.

We will use a shortened notation $g_{i}=g(t,\omega_{i},s_{i})$, so $Q=q\cdot g$. Note also $\ddot{g}=-\lambda\cdot g$. We can now write the impact conditions from Eq(1) in a matrix form $Ab=0$, where

$g=g(\tau)$, $g^{\prime}=g^{\prime}(-\tau^{\prime})$ and $b=[q;q^{\prime};F^{0}_{N}]$. $A$ is a $(2N+1)\times 2N$ matrix. For a nontrivial solution $b$ to exist, $A$’s rank must be at most $2N-1$. To evaluate $A$’s rank, it is convenient to bring it to a block upper triangular form by applying rank-preserving transformations. Using $m^{-1}kX=\lambda^{\mathsf{T}}\cdot X$, $cX^{\mathsf{T}}mX=I$, $\eta^{\mathsf{T}}=X_{N}\cdot X_{N}$, $\ddot{g}=-\lambda\cdot g$ and Eq(8), after straightforward manipulations, one finds (see B for details):

where

The condition $\textrm{rank}(A)<2N$ is equivalent to $\textrm{rank}(B)<N$, which under the assumption of non-singular $\bar{U}$ can be written as $\det{B_{(N)}}=\det{B_{(N+1)}}=0$ in the form:

We will call the above equations, written in this or an equivalent form, the impact equations. We would like to stress that the impact equations only involve $\tau^{p}$ and $\lambda^{p}$, thus validating our claim that the collisionless solution is determined by the spectra $\lambda^{p}$ alone.

As was explained in Section 2, the collisionless solution is invariant to $TS^{p}$ associated with the symmetry point $P^{p}$, where $S^{p}$ is a Euclidean metric preserving transformation, i.e. a reflection or rotation. For a non-degenerate spectrum each normal mode must respect $TS$ symmetry, leaving only two possibilities at $P$: a) $q_{i}$ is unchanged by $S$ and $g_{i}(t)$ is symmetric ($|s_{i}|=1$), b) $q_{i}$ is flipped by $S$ and $g_{i}(t)$ is antisymmetric ($s_{i}=0$). Without loss of generality, we will consider non-negative $s_{i}$. In that case, $g_{i}(t)/\dot{g}_{i}(t)=\sigma_{i}\tan^{\sigma_{i}}(\omega_{i}t)/\omega_{i}$, where $\sigma_{i}=1-2s_{i}$. Therefore,

Eq.(15) is easily solvable for $N=2$. For $G$ one finds:

Using Eq.(16), we can write the impact equations as

Since the symmetry points are separated by an impact, we are interested in the solution with $tt^{\prime}<0$, (correspondingly, $\tau\tau^{\prime}>0$). Note also, for $\omega^{2}<0$

We found that the impact equations have a non-trivial solution (i.e. a solution with finite impact times) only for

see C for details.

It is convenient to express the impact equations using dimensionless notations $o_{i}^{p}=|\omega_{i}^{p}|\tau^{p}$, which were called impact phases in [5]. For example, for $\lambda_{i}>0$, Eq.(18) reads in terms of the impact phases

where $\mu=\tau/\tau^{\prime}$. If $\omega^{2}<0$, then $\sigma\tan^{\sigma}(o)o$ should be replaced with $-\tanh^{\sigma}(o)o$, according to Eq.(19) .

Several collisionless models with $N=2$ have been considered in the literature. We show below how their corresponding impact equations are captured by Eq.(18). Their solutions can be conveniently expressed in terms of the positive roots of the equation

defined such that $y_{n}(a,b)\in[(n-1)\pi,n\pi)$. We also define $\beta_{n}(\rho)=y_{n}(-\rho,\rho)$, $\gamma_{n}(\rho)=y_{n}(1/\rho,\rho)$ and $\alpha_{n}=\gamma_{n}(0^{+})$. We will list the solutions in terms of the impact phases $o_{2}$ and $o^{\prime}_{1}$, representing the impact times $\tau=o_{2}/\omega_{2}$ and $\tau^{\prime}=o^{\prime}_{1}/\omega^{\prime}_{1}$.

For $N=2$ we obtained analytical solutions (expressed via univariate functions $\beta_{n}$ and $\gamma_{n}$) of the impact equations for various symmetry point configurations $\sigma^{p}$, and proved the solution existence condition Eq.(29). While for $N>2$ a closed-form solution might not be available, we can still analyze the impact equations and provide compelling arguments that the same solution existence condition remains valid for all $N$. In addition, in this section, we present a number of analytical results obtained for vanishing $\lambda^{\prime}_{N-1}$.

We define

and hence

One can observe that: a) for $\lambda^{\prime}_{N-1}=0$ and $G^{N-1}=0$, we have $\textrm{rank}(B)<N$, b) all quantities entering $B$, apart from $w^{\prime}$, are analytic in $\lambda^{\prime}$, while $w^{\prime}$ is analytic in $\omega^{\prime}$. It then follows that $w^{\prime}_{N-1}$ is non-singular along the solution of the impact equations for $\lambda^{\prime}_{N-1}\to 0$, i.e.

where $c_{0}$ is some constant. With this insight, it is straightforward to derive equations on $\tau^{p}$ for $\lambda^{\prime}_{N-1}\to 0^{+}$ by differentiating $\det{B_{(N)}}$ and $\det{B_{(N+1)}}$, as we explain below.

First, let us show that if $\lambda^{\prime}_{N-1}=0$ and $G^{N-1}=0$, then $\textrm{rank}(B)<N$. Let $\tau^{p}(\lambda^{\prime}_{N-1})$ represent the dependence of the impact equations’ solution on $\lambda^{\prime}_{N-1}$, with the rest of the spectra being fixed. For a quantity $A$ depending on $\tau^{p}$ and $\lambda^{\prime}_{N-1}$, let $A(\lambda^{\prime}_{N-1})$ be a shorthand notation for $A\left(\tau(\lambda^{\prime}_{N-1}),\tau^{\prime}(\lambda^{\prime}_{N-1}),\lambda^{\prime}_{N-1}\right)$. We define $\lambda^{\prime}_{N-1}=0$ limit as:

which we denote with a left superscript $0$. If we assume ${\,}^{0}G^{N-1}=0$ then ${\,}^{0}U^{N-1}={\,}^{0}\!M^{N-1}=1/\lambda$. Therefore, ${\,}^{0}\!B^{N-1}={\,}^{0}\!B^{N}$, and hence $\textrm{rank}({\,}^{0}\!B)<N$. From Eqs.(30,32) it follows (for $c_{0}>0$, which remains to be shown) that $\tau^{\prime}\to+\infty$ for $\lambda^{\prime}_{N-1}\to 0^{+}$, and hence ${\,}^{0}G$ and ${\,}^{0}U$ are independent of $\tau^{\prime}$. We also define a differential operator

Expanding $\det{B_{(N)}}=\det{B_{(N+1)}}=0$ around $\lambda_{N-1}^{\prime}=0$, one concludes that the necessary condition for $\textrm{rank}(B)<N$ becomes, to first order,

Let a matrix $\tilde{B}$ have the following properties: a) $\tilde{B}$ is $N\times N$ matrix, b) ${\,}^{0}\!\tilde{B}^{N-1}={\,}^{0}\!\tilde{B}^{N}$, c) $\bar{\tilde{B}}=[\bar{U},1/\bar{\lambda}]$. In the rest of the chapter, to simplify notations, we will be omitting the left superscript wherever the limit of Eq.(33) is clearly implied. One can prove (see D for details) the following formula:

From Eq.(36) we find for $\tilde{B}=B_{(N)}$

and for $\tilde{B}=B_{(N+1)}$

We define

Note that $K$ and $\tilde{K}$ do not depend on $\tau^{\prime}$, (since ${\,}^{0}U$ only depends on $\tau$) . From Eqs.(37,38,39) we see that the conditions in Eq.(35) can be written as $(K+\tilde{K})[1;1/c_{0}]=0$, from where we have the following (essentially decoupled) equations on $\tau$ and $\tau^{\prime}$:

For $N=2$ these equations further simplify to

which is in agreement with the results of Section 4.

For a general $N$, it is also straightforward to solve for $\tau$ in the asymptotic case of large $\tau$, when the number of oscillations in the unconstrained phase is large. See D for details.

We have shown that for $\lambda^{\prime}_{N-1}\to 0$, Eq.(32) holds for any $N\geq 2$. If $c_{0}>0$, it follows that $o^{\prime}_{N-1}$ is constant in the limit $\lambda^{\prime}_{N-1}\to 0^{+}$ (see Eq.(63) for an explicit formula) and hence $\tau^{\prime}\propto 1/\omega^{\prime}_{N-1}$, while no $\tau^{\prime}>0$ solution exists for $\lambda^{\prime}_{N-1}\to 0^{-}$. For $N=2$, $c_{0}=-w_{1}=\tanh^{\sigma_{1}}(o_{1})\nu_{1}>0$. For $N>2$ we empirically verify the condition $c_{0}>0$ by randomly sampling the spectral values. While not a rigorous proof, it strongly suggests the correctness of the general collisionless solution existence condition Eq.(29).

We found it helpful to visualize the solution of the impact equations using a contour-plot in $(o_{N},o^{\prime}_{N-1})$ coordinates (equivalent to $(\tau,\tau^{\prime})$ coordinates, up to a rescaling). We could plot a function that turns zero whenever either of the impact equations is satisfied, such as $\det{B_{(N)}}\det{B_{(N+1)}}$. However, this function is discontinuous (where $\dot{g}_{i}=0$ or $g^{\prime}_{i}=0$). Instead, we will use

which is continuous, and for which $\textrm{rank}(\check{B})=\textrm{rank}(B)$ (away from the discontinuities of $B$). In Fig.(2) we plot the zero contour of

Different curved lines correspond to the solutions of $\det{B_{(N)}}=0$ and $\det{B_{(N+1)}}=0$, while their intersections represent the solutions of the impact equations. Whether such a solution is physically realizable needs to be further verified by ensuring that no unwanted ground penetration happens in the unconstrained phase and no ground contact loss happens in the constrained phase. In the figure we consider a near critical configuration of a randomly sampled (and shifted) spectrum, so that $0<\lambda^{\prime}_{N-1}\ll-\lambda_{N-1},\lambda_{N}$ for $N=4$. In that regime the picture is qualitatively similar for different values of $N$, as also follows from our analysis of the critical region (in this section and D). We also depicted the analytical asymptotic large-$\tau$ solution presented in D, see Eq.(63). The asymptotic solution forms a square grid with the step size $\pi$. We only depict the bottom row of the grid, relying on the intuition from $N=2$ case, for which one can verify that only the bottom row is physically realizable.

So far, we have developed a general approach for efficiently finding collisionless solutions in complex models. We have reduced the problem to a nonlinear equation in just two variables (regardless of the model complexity) – the impact times. We also proposed a solution existence condition and showed how $N=2$ collisionless solutions, considered in the literature, easily follow from our formulation. In this section we illustrate our approach by considering a rocking motion of a planar 3-DOF biped (i.e. $N=3$), which to our knowledge has not been previously studied.

The model consists of three parts (see Fig.(3)): 1) a pair of legs, rigidly affixed together at an angle $2\theta$, 2) a standing torso, attached to the legs at the hip, and 3) a hanging arm, attached at the top of the torso. For simplicity, all the links have the same length $l$. Again, we consider a harmonic approximation, which becomes exact in the limit $\theta\to 0$.

It has been suggested in the literature to use torsion springs to prop up a standing torso in the context of collisionless walking models [14, 4, 5]. In ref [5] it was also hypothesized that such springs are not strictly necessary for keeping the torso pointed up, as long as the torso is endowed with a hanging arm. An explicit collisionless solution, that we obtain below, can be viewed as a definitive proof of that hypothesis. In fact, this result holds for a torso containing any number of links stacked up on top of each other. This directly follows from our solution existence condition (without the need of an explicit solution), by noticing that a single hanging link implies $\lambda^{\prime}_{N-1}>0$.

The angle conventions and masses of the model are shown Fig.(3). For the masses we use unitalicized $\mathrm{m}_{i}$ to distinguish them from the mass matrix $m$. The total mass of the model is $\mathrm{m}=\mathrm{m}_{f}+\mathrm{m}_{1}+\mathrm{m}_{2}+\mathrm{m}_{3}$, where $\mathrm{m}_{f}=2\mathrm{m}_{0}$ is the mass of the “feet”. The coordinates $x=[x_{1};x_{2};x_{3}]$ are the angles of the arm, torso and the stance leg, measured relative to the static equilibrium configuration $x^{eq}$, which (when measured relative to the vertical axis) is $x^{eq}=[\pi;0;\xi]$, where $\xi=\theta\mathrm{m}_{f}/(\mathrm{m}-\mathrm{m}_{f})$. Note that $F_{N}^{0}=\theta\mathrm{m}gl$ and $x^{0}=[0;0;x_{N}^{0}]$, where $-x_{N}^{0}=\theta+\xi=\theta\mathrm{m}/(\mathrm{m}-\mathrm{m}_{f})$.

Similar to $N=2$ rocking motion, considered in Section 4.3, the symmetries are: $\sigma=[-1;-1;-1]$ and $\sigma^{\prime}=[1;1]$.

The mass $m$ and spring constant $k$ matrices are:

and

In the constrained phase $m^{\prime}=m_{(NN)}$ and $k^{\prime}=k_{(NN)}$.

Let us summarize the formal steps of finding a collisionless solution. First, solve for the eigensystem of $(m^{p})^{-1}k^{p}$ to find $\lambda^{p}$ and $X$ (see the definition of $X$ near Eq.(4)). Construct $M$ (see Eq.(7)) and find $X^{\prime}$ (see Eq.(8)). Construct $\check{B}$ (see Eq.(42) for $\check{B}$ and Eqs.(10,11) for $\eta$ and $g^{p}$ to that end). Plot the zero contour of $\phi(o_{N},o^{\prime}_{N-1})$ (see Eq.(43)). Use the (approximate) locations of curve intersections in the plot as an initial guess of the solution of $\det{\check{B}_{(N)}}=\det{\check{B}_{(N+1)}}=0$ to solve that system numerically to find $\tau^{p}$. Given the knowledge of $\tau^{p}$, $g^{p}$ and $X^{p}$, construct $A$ (see Eq.(12)). (Optionally, as a correctness check, verify that $\textrm{rank}(A)<2N$). Finally, use the top left corner $(2N-1)\times(2N-1)$ submatrix of $A$ to find $q^{p}$ as²²2$q^{p}$ can be expressed in terms of $\tilde{A}$ and $x^{0}$ without resorting to $X^{p}$. Although interesting theoretically, this is irrelevant numerically.

The collisionless trajectory is then given by $x=Xq\cdot g$ in the unconstrained phase and $x=X^{\prime}q^{\prime}\cdot g^{\prime}+x^{0}$ in the constrained phase.

The described procedure is computationally trivial, taking a few seconds of running an Octave script [15], with almost all of the time spent on the contour plot, needed for a reasonably good initial guess of the impact equation solution.

The rest of the results in this section are presented for $g=l=\mathrm{m}_{i}|_{i=0,...,3}=1$. Then $\xi=2\theta/3$, $F_{N}^{0}=5\theta$ and $x_{N}^{0}=-5\theta/3$. The zero-contour plot of $\phi(o_{N},o^{\prime}_{N-1})$ is shown in Fig.(4). Its qualitative similarity to Fig.(2) is not surprising, as in both cases each phase has exactly one positive eigenvalue $\lambda^{p}_{N}$, due to a hanging arm in this model. Again, the impact equation solutions appear to form a square grid asymptotically for large $\tau^{p}$. In the grid, columns and rows correspond to different number of oscillations in the unconstrained and constrained phases, respectively. Again, only the bottom row is physically realizable, as the impact equation solutions in the higher rows have stretches of time (in the constrained phase) with an attractive ground reaction force. We investigate a solution with the smallest number of oscillations, corresponding to the left-most intersection in the bottom row, marked with a circle in Fig.(4).

Solving for the impact times, as described above, we find $t_{imp}=\tau=3.0795$ and $t^{\prime}_{imp}=-\tau^{\prime}=-0.77785$. See E for numerical values of other quantities. The complete solution $x^{p}(t)$ in the units of $\theta$ is shown in Fig.(5), (also showing $\dot{x}^{p}(t)$ and $\ddot{x}^{p}(t)$). We see that, as was described in Section 2, $x$, $\dot{x}$ and $\ddot{x}$ are all continuous at the impact, with $\dot{x}_{N}$ and $\ddot{x}_{N}$ vanishing, and $\ddot{x}$ displaying a cusp. In Fig.(6), the collisionless solution at $P$, impact, and $P^{\prime}$ is depicted to scale for $\theta=0.1$.

In this work we have extended analytical analysis of collisionless motion to $N$-DOF models with one-dimensional impacts. The analysis revealed that the possibility of expressing the impact equations in terms of the impact phases is a general property, not limited to small $N$. The relative simplicity of this formulation enabled the investigation of the solution in the critical region, demarcated by the proposed solution existence condition.

With this development, a host of new models (e.g. multi-DOF extensions of $N=2$ models discussed in Section 4) can be further investigated. However, to be able to apply this formalism to general two-dimensional (e.g. the model of [4]) and three-dimensional bipedal walkers, it must be further extended to cover two- and three-dimensional impacts, respectively. Our preliminary consideration suggests that in addition to relying on the constrained and unconstrained spectra, it will also involve partially-constrained spectra – one per each dimension of the impact.

We would like to point out how infinite friction could be used to facilitate solving for collisionless motion with higher-dimensional impacts. As was explained in [5], a collisionless solution for $(n+1)$-dimensional impact requires adjustment of $n$ model parameters. For example, for a one-dimensional impact $n=0$ (as in this paper), so a collisionless solution can be found for arbitrary model parameters, as long as the solution existence condition is met. Assume now that $n>0$ and the acceleration vanishing (at the point of contact, e.g. foot-ground contact) is only enforced for the component orthogonal to the ground. In that case, the foot will slip right after touching the ground, and the amount of slippage will be inversely proportional to the friction coefficient $\mu_{f}$. Clearly, the slippage is eliminated for $\mu_{f}\to\infty$. Consequently, in the presence of infinite friction, there is a collisionless solution for arbitrary model parameters even for $n>0$. In that case, collisionless solutions are easily found using the contour-plot method described in Section 5. A suitable solution (e.g. with minimal number of torso oscillations) can then be continuously evolved numerically [5] back to finite $\mu_{f}$, with the desired topological properties of the solution preserved.

We believe that the symmetries imposed in this work are not essential for collisionless motion; they were used because of the simplifications they bring. Finding a collisionless solution without reliance on the above symmetries would be a step toward considering more faithful anthropomorphic and animal models.

While more research remains to be done, we believe that the design of complex models with desired anthropomorphic (as well as non-anthropomorphic) form factors, capable of perfectly energy-conserving mechanical motion (i.e. with zero COT), is within a reach. This paper is a contribution in that direction.

We thank Georges Harik for many useful discussions.