A minimalistic stochastic dynamics model of cluttered obstacle traversal

by   Bokun Zheng, et al.
Johns Hopkins University

Robots are still poor at traversing cluttered large obstacles required for important applications like search and rescue. By contrast, animals are excellent at doing so, often using direct physical interaction with obstacles rather than avoiding them. Here, towards understanding the dynamics of cluttered obstacle traversal, we developed a minimalistic stochastic dynamics simulation inspired by our recent study of insects traversing grass-like beams. The 2-D model system consists of a forward self-propelled circular locomotor translating on a frictionless level plane with a lateral random force and interacting with two adjacent horizontal beams that form a gate. We found that traversal probability increases monotonically with propulsive force, but first increases then decreases with random force magnitude. For asymmetric beams with different stiffness, traversal is more likely towards the side of the less stiff beam. These observations are in accord with those expected from a potential energy landscape approach. Furthermore, we extended the single gate in a lattice configuration to form a large cluttered obstacle field. A Markov chain Monte Carlo method was applied to predict traversal in the large field, using the input-output probability map obtained from single gate simulations. This method achieved high accuracy in predicting the statistical distribution of the final location of the body within the obstacle field, while saving computation time by a factor of 10^5.


page 1

page 6

page 8


Environmental force sensing enables robots to traverse cluttered obstacles with interaction

Many applications require robots to move through terrain with large obst...

Shape-induced obstacle attraction and repulsion during dynamic locomotion

Robots still struggle to dynamically traverse complex 3-D terrain with m...

Geometry and Dynamics for Markov Chain Monte Carlo

Markov Chain Monte Carlo methods have revolutionised mathematical comput...

Deep reinforcement learning with a particle dynamics environment applied to emergency evacuation of a room with obstacles

A very successful model for simulating emergency evacuation is the socia...

Lambda-Field: A Continuous Counterpart Of The Bayesian Occupancy Grid For Risk Assessment And Safe Navigation

In the context of autonomous robots, one of the most important tasks is ...

Bayesian Approach to Inverse Time-harmonic Acoustic Scattering from Sound-soft Obstacles with Phaseless Data

This paper concerns the Bayesian approach to inverse acoustic scattering...

A Tractable Analysis of the Blind-spot Probability in Localization Networks under Correlated Blocking

In localization applications, the line-of-sight between anchors and targ...

I Introduction

Many important applications require robots to traverse cluttered obstacles such as search and rescue in rubble, environmental monitoring on the forest floor, and extraterrestrial exploration through Martian rocks. The dominant approach of robotic locomotion in complex environments is to avoid obstacles [4, 3, 17], which could be challenged in such cluttered terrain. By contrast, insects such as the discoid cockroach are excellent at traversing cluttered large obstacles, often using their bodies and appendages to make physical interaction with obstacles [11, 20, 10, 6, 7, 9] (for a review, see [16]). Understanding how the physical interaction between a self-propelled locomotor and obstacles leads to or hinders traversal can help inform robot locomotion beyond avoiding obstacles and increase the terrain that they can access.

Recent animal and robot studies from our lab established a potential energy landscape approach for quantifying locomotor-obstacle interaction, which provided insights into the mechanical principles of how different modes of locomotion emerges from the interaction, which can be controlled to achieve traversal [16, 15]. However, there still remains a lack of understanding of the dynamics of cluttered obstacle traversal, which is not modeled by the potential energy landscape approach. Here, we take the next step towards this by developing a dynamic model of locomotor-obstacle interaction using a minimalistic model system and studying how the properties of locomotor self-propulsion and obstacle affect traversal.

Our study is inspired by and builds upon findings from our previous work [15]. The potential energy landscape approach revealed that, when interacting with large obstacles, a self-propelled locomotor has highly dissipated dynamics [16] and is strongly attracted to basins of a potential energy landscape resulting from locomotor-terrain interaction [16, 15]. To escape from attraction to certain basins that lead to unfavorable modes (e.g., become trapped against the obstacles [9, 15]), the locomotor can use both active sensing and control and passive kinetic energy fluctuation from self-propelled oscillation [15]. As a first step towards understanding traversal dynamics, we focus on stochastic dynamics from self-propelled oscillation and do not consider active sensing and control. In addition to adding dynamics, our work further explores the interplay of self-propulsion and random oscillation and the effect of asymmetry in terrain physics common in nature [12, 19] and how locomotor-obstacle interaction leads to traversal over large spatiotemporal scales, which the previous work did not address.

Our minimalistic model system consists of a circular locomotor body self-propelled forward with only two translational degrees of freedom moving on a level plane interacting with a pair of grass-like beam obstacles that rotate within the plane. A lateral random force is added to introduce body oscillation, leading to stochastic dynamics. We choose to use a random force considering that animals often have random motion during locomotion

[5, 2, 21, 13]. In addition, we introduce asymmetry in the environment by allowing the two beam obstacles to have different torsional stiffness. Moreover, we used simulation results with systematically varied propulsive and random forces and asymmetric beam stiffness to test how well the potential energy landscape can qualitatively predict traversal dynamics without solving equations of motion. Furthermore, we extended the single locomotor-obstacle interaction model to a larger terrain with multiple obstacles in a lattice configuration test how well it can statistically predict how the body moves across a complex terrain on a larger scale.

We introduce how the model was defined and the dynamics modeled (Sec. II-A, II-B), how the potential energy landscape was calculated (Sec. II-C), and how we generated a multi-obstacle field (II-D). We then describe results using the stochastic dynamics simulation to study traversal of a pair of beams and testing how well the landscape model informed traversal (Sec. III) and how well the traversal of a multi-obstacle field can be predicted by a Markov chain Monte Carlo method (Sec. III-C). Finally, we summarize our findings and their implications (Sec. IV).

Ii Methods

Ii-a Model Definition

Our simplistic 2-D model of body-obstacle interaction (Fig.1) consists of three rigid bodies on a horizontal plane: a circular disc representing the locomotor body (mass = 1 kg, radius = 10 m), and two plates with viscoelastic revolute joints at their far ends (mass = 0.1 kg, length

= 25 m, moment of inertia about joint

= 20.83 kgm). Hereafter, we refer to them as beams. The beams align with each other when unloaded, forming a closed “gate”. For simplicity, we assumed that the body surface is frictionless. Thus, body-beam contact results only in forces normal to the body surface, which generates no torque. Thus, body orientation remains unchanged. In addition, we assumed that 2-dimensional body-beam collisions are inelastic with partial energy loss. We used coefficient of restitution to measure the degree of elasticity, defined as the ratio of final to initial relative normal velocities between the body and beams during each collision. Here we tuned by observing different trials and chose a value of 0.8 to conserve a large portion of mechanical energy so that traversal still occurs even after a large number of collisions. Finally, for simplicity, we assumed that there is no friction against the ground (although it can be added–see Sec. III-C, which should only quantitatively but not qualitatively change our results).

Fig. 1: Schematic of the 2-D model layout A screenshot of the simulation at initial state. The circle represents the body, the blue/red line represents the left/right beam, the orange arrow represents the propulsive force, and the gray arrow represents the random force that oscillates laterally.

Besides forward propulsion, animals with a sprawled leg posture like cockroaches can generate substantial lateral forces during locomotion [18]. In addition, animals can randomly change their movement direction during locomotion [13]. Considering these and for simplicity, self-propulsion of the locomotor body was modeled by a constant forward propulsive force and an oscillating random lateral force (Fig. 1

red and green arrows) with a standard Gaussian distribution:


where is a constant characterizing the magnitudes of the random forces and

, a normal distribution with a mean of 0 and a standard deviation of 1.

Because the body does not rotate, forward and lateral forces are always in the - and - directions, respectively. Applying Newton’s second law in these two directions, the equations of motion of the body are:


where is body mass, is body acceleration, and is the resistive force from the beam, with for the left and right beams, respectively.

Applying Newton’s second law to the beam:


where is beam moment of inertia about its joint, are the orientation, angular velocity, and angular acceleration of the beam, and and are the torsional stiffness and damping coefficient of the beam joint.

Ii-B Interaction Dynamics

We used the Euler method to numerically integrate forward in time to calculate the dynamics of the system. A time step of 0.004 s was chosen to achieve good numerical accuracy while maintaining computation efficiency. The run-time of each trial varied from 5 to 20 min on a lab workstation. The key part of dynamics is to determine the interaction between the body and beams. We developed two complementary methods to model interaction dynamics for two different scenarios: a collision method and a constraint method. The collision method models repeated collisions between body and beams after initial contact, which occur due to the competing forward propulsive force and backward beam resistive forces. However, as these collisions dissipate energy, the duration of single collision decreases and eventually becomes smaller than the time step, leading to a significant increase in numerical error [1]. In this case, we assumed that the interacting body and beam have no relative motion in the direction normal to the contact point but can move relatively in the tangential direction and developed the constraint method to describe interaction dynamics.

Ii-B1 Collision method

When calculating collision occurring over an infinitesimal time, we considered the body and beam as a system with the center of rotation at the beam joint. The system is subject to a finite external beam joint torque and external forces . After multiplying the finite values with an infinitesimal time, the external momentum impulse is negligible. In addition, there is a reaction force acting on the fixed beam joint, which has a zero moment arm. Therefore, the angular momentum of the system about the instantaneous center of rotation is conserved during a collision.

The calculation of collisions depends on the geometric relationship between body and beam. During traversal, there are two different contact cases: tangential contact and point contact.

Fig. 2: Schematics of the tangential contact between the body and left beam before and after a collision.

In the tangential contact case (Fig. 2), the beam is tangential to the body surface at the point of contact. The normal direction of the body surface at the point of contact is perpendicular to the beam. Applying conservation of angular momentum to the body-beam system:


In addition, applying the definition of CoR to the tangential contact case:


where is the distance from beam joint to contact point, and are beam angular velocities before and after collision, and are local velocities of the beam at the contact point before and after collision, and are body velocities before and after collision, , , , and are their projections in the direction normal and tangential to the beam. Note that = because there is no interaction force in the tangential direction from our frictionless body assumption. Using Eqns. 4, 5, we can solve the dynamics of the collision in the tangential contact case.

Fig. 3: Schematics of the point contact case before and after a collision.

In the point contact case (Fig. 3), the free end of the beam contacts the body surface. Because the normal direction of the body surface at the contact point is no longer perpendicular to the beam, it is more complex to formulate a expression in scalar form. Considering that its collision is small and only happens in the short duration before detachment, we assumed the collision to be perfectly elastic in this case for simplicity. Thus, mechanical energy is assumed to be conserved at the collision time step:


Also the system complies with the conservation of angular momentum before and after collision:


where is the change of body momentum, and are body velocities before and after the collision, is the change of angular momentum of the beam before and after the collision, is the interaction force acting on the body from the beam, is the distance from the beam joint to contact point, and is the change of angular momentum of the system, which is zero at collision. The velocities after the collision can be solved from those before using Eqns. 6-9.

Ii-B2 Constraint method

The constraint method is based on Newton’s second law and the geometric constraints that the body and beam have the same normal velocity at the contact point. It also differs between the two contact cases.

In the tangential contact case, the velocity constraint is:


Taking time derivative of both sides of Eqn. 10, we have:


where is the center of mass position measured from the beam joint, is body velocity, and is the body acceleration.

In the point contact case, the velocity constraint is:


Taking time derivative of both sides of Eqn. 12, we have:


where is beam length.

In the simulation, the collision and constraint methods are integrated following the workflow shown in Fig. 4. The simulation first runs the collision method in loops. As the change of the body’s momentum during each collision decreases below a small threshold kgm/s, i.e., collisions become small enough and exchange little momentum, the simulation switches to the constraint method in loops. During the switch, beam angular velocities are updated to satisfy the initial velocity constraint of the constraint method.

Fig. 4: A flowchart of the simulation program.

Ii-C Potential Energy Landscape

The potential energy landscape is a modeling approach to model the conservative forces during locomotor-terrain interaction over relevant degrees of freedom. Without knowledge of non-conservative forces that are often difficult to measure, it provides a useful approach for understanding how the system may or may not move on the landscape, as long as potential energy landscape plays a major role in shaping dynamics [16]. In our system, potential energy arises from elastic torsional joints of the beams:


where is the total beam potential energy, is the elastic stiffness of the L/R beam, and is the deflection angle of the beam, with for the left and right beams, respectively.

The landscape is calculated as a 3-D surface over a mesh grid of m. At each grid point, we assumed that beams are always deflected forward and contact the body if the body is within the range of beam radius. During traversal, the body pushes forward and deflects the beams, resulting in two potential energy barriers on each side (Fig. 5AB). The two barriers overlap in the middle area, where the body can interact with both beams. Fig. 5C is an example of asymmetric landscape with different and .

Fig. 5: Potential energy landscape for the beam obstacle traversal. (A) Landscape with only the elastic energy. = 400 Nm/rad. (B) A top view of the landscape in A. (C) An asymmetric landscape with = 500 Nm/rad and = 250 Nm/rad. (D) Landscape with sum of elastic energy and propulsive conservative energy. = 400 Nm/rad.

Because we assumed a constant forward , the potential energy landscape can also include conservative potential energy from . The summed potential energy is:


where is the potential energy from . The zero level of was defined at the anterior boundary = 60 m.

A basin that spans across the direction in front of the beams (around = 20 m) emerges after including the energy of of a modest magnitude in the potential energy landscape (Fig. 5D). This helps us understand how the landscape is tilted by the propulsive force. Within insufficient , the locomotor body would be trapped in this basin and not traverse. Adding kinetic energy fluctuation induced by the random force can help the body escape the horizontal basin. With sufficient , the landscape is so heavily tilted that the basin disappears, and the locomotor should always traverse.

The evolution of the potential energy of the system in the simulation traversal trials can be plotted over the landscape to visualize how the body-beam interaction influences body motion. Typically, the trajectory of the system potential energy is above the landscape surface due to the beam inertia. However, sometimes the trajectory penetrates the landscape surface. An example is shown in Fig. 6B. In this case, the left beam lost contact with the body and returned its initial orientation (Fig. 6

A). Because the body is still within the radius of the beam, our landscape calculation that always assumes beams being deflected forward over-estimates the potential energy. If we only consider the landscape from the right beam that is deflected, this artifact is removed (Fig.

6C). See multimedia material part 3 for an example video.

Fig. 6: An example of trajectory penetrating the landscape surface. Points with show the mapping between 2-D position and landscape. (A) The screenshot in simulation. (B) The trajectory plotted over the normal landscape. (C) The trajectory plotted over the landscape without the inactive barrier.

Ii-D Multi-Obstacle Field Interaction

Using the body-beam interaction model above with a single pair of beams (a gate), we can create a large cluttered obstacle field by composing multiple pairs of beams (gates) in a lattice arrangement. Fig. 7 shows a 9 9 gate obstacle field used for further simulation. Except for the first gate region which uses a set of initial system states (), the system state input of each gate region comes from the output of the previous gate region. An index was used to track which gate region is being activated. For this 9 9 gait obstacle field, and . Finally, the main program outputs a trajectory restricted in a single gate region, which is a series of local positions , where m, m. Later in the visualization process, according to the gate index , are converted to the global position data , where m, m. See multimedia material part 4 for a simulation example.

Fig. 7: A 9 9 lattice obstacle field. Green dashed line show boundary of the region of each pairs of beams. show beam joints. Blue and red line segments show left and right beams of each pair.

The pseudo code:

while  and  do
     Single gate simulation,
     if Reached left boundary, then
     else if Reached right boundary, then
     else if Reached top boundary, then
     end if
end while
function mat2avi()
     Convert local position data,
     return video.avi
end function
Algorithm Gate lattice simulation

All the modeling calculations and simulations were performed in MATLAB R2019b.

Iii Results

Iii-a Traversing symmetric beams

We first studied the probability of traversing a pair of symmetric beam obstacles under self-propulsive and random forces. The body starts at an initial forward velocity , accelerates for a short distance before making initial contact with beams. The criterion of traversing successfully is that the body can reach = 65 m, where the body is out of the range of beams (multimedia material, part 1). Failure resulted from: (1) the body gets stuck by the beams (multimedia material, part 2), (2) the body deviates from the desired track and exits from the side boundaries ( -25 m or 25 m), and (3) the body did not reach = 65 m within maximal iterations.

From the potential energy landscape (Fig. 5A), the locomotor body needs to have sufficient kinetic energy to overcome a potential energy barrier to traverse [15]. Because higher propulsive force can tilt the landscape, facilitating to overcome barriers by lowering them, while the random force adds more kinetic energy fluctuation in the system. Thus, we hypothesized that the probability of traversing a pair of symmetric beam obstacles increases with propulsive force and random force.

To test this hypothesis, we varied propulsive force and the magnitude of random force using system parameters in Table I and ran 100 trials for each given () to obtain traversal probability (Fig. 8). For a given , traversal probability increased with propulsive force. However, for a given , traversal probability did not always increase with random force. A random force of up to 20 N increased traversal probability. However, as increased beyond 20 N, traversal probability decreased. Examination of simulation videos revealed that this was because more trials exited the side boundaries of the gate region before arriving at the beams. Thus, a sufficient but not excessive lateral random force helps traverse symmetric beam obstacles.

Parameter Value
Propulsive force 4,5,6,7,8,9 N
Random force magnitude 0,10,20,30,40 N
Beam stiffness 400 Nm/rad
Beam damping 50 Ns/rad
Initial forward velocity 0 m/s
Coefficient of Restitution 0.8
Oscillator frequency 50 Hz
Acceleration distance 20 m
Max iterations 3000
TABLE I: Symmetric test configuration
Fig. 8: Probability of traversing a pair of symmetric beams as a function of propulsive force and random force magnitude.

Iii-B Traversing asymmetric beams

Next, we studied the probability of traversing a pair of asymmetric beams with different stiffness. The potential energy landscape becomes asymmetric in this case, with a barrier lower on one side where the beam has a smaller stiffness (Fig. 5C). We hypothesized that the probability of traversing on the side of lower stiffness increases with the level of asymmetry.

To test this hypothesis, we varied the left beam stiffness while keeping the right beam stiffness constant and varied and systematically (Table II). Unlisted parameters are the same as those in Table I. Here we selected three groups with varying stiffness ratios as examples to illustrate the typical results (Fig. 9). Only successful traversal trajectories are plotted over the landscape. In the symmetric control group with both beam stiffness being 500 Nm/rad (Fig. 9A), traversal was rare with similar probability on both sides (L:13%, R:8%). This is because the propulsive force was barely sufficient to overcome the potential energy barriers. When the left beam stiffness was reduced from 500 to 300 Nm/rad, traversal was frequent (52%), with the body pushing across the left beam (Fig. 9B). As the left beam stiffness further reduced to 100 Nm/rad, traversal was almost guaranteed by pushing across the left beam (98%, Fig. 9C).

Parameter Value
Propulsive force 7,8N
Random force magnitude 10,20,30,40,50,60 N
Left beam stiffness 100:50:500 Nm/rad
Right beam stiffness 500 Nm/rad
TABLE II: Asymmetric test configuration
Fig. 9: Successful trajectories on the landscape in asymmetric test. Self-propulsive force = 7 N, random force magnitude = 10 N, Nm/rad. (A) Nm/rad. (B) Nm/rad. (C) Nm/rad.

For all the asymmetric tests varying and

, we classified the trajectories into “left” or “right” types, reflecting traversal tending to which barrier side, by comparing the maximal deflected angle of the left or right beam. Then, we can quantitatively evaluate the asymmetry of trajectories by calculating the ratio of the number of trials for each type to the total number of successful traversal trials (Fig.

10). For N, as the left beam stiffness decreases, traversal becomes less frequent on the right side. However, for N, this trend is still existed, but less monotonically due to the high stochasticity. In addition, traversal probability increased with (Fig. 10A vs. B). These results well supported our hypothesis and were consistent with observations in the traversal of symmetric beams.

Fig. 10: Ratio of right trials to all successful trials as a function of left beam stiffness and random force magnitude. A darker block in the matrix indicates that there are fewer right trials than the left ones, and traversal becomes more asymmetric towards the left. A brighter block indicates the number of left and right trials are more equal.

Iii-C Multi-Obstacle Traversal

On the larger terrain, trajectories become more complicated and lack of patterns. Here we focused on how the locomotor body transits between adjacent gate regions. As a first step, we assumed that all the gates are symmetric with the same stiffness (400 Nm/s). Each gate is considered as a single system with an input and an output when crossing region boundaries. Here we discretized the boundary of a gate region into 6 segments: middle traverse (MT), left traverse (LT), right traverse (RT), left deflect (LD), right deflect (RD) and enter (EN) (Fig. 11). For the traversal problem, EN is only for input and MT is only for output, whereas the other four can be either input or output. Another possible output is that the body is trapped inside a gate region in between the two beams. Motion in each gate can be treated individually as a stochastic process independent of the history before entering. And we assume for fixed model parameters, there is a constant probability correspondence between input and output states. Thus, the Markov Chain [14] can be applied to describe the stochastic transition between gate regions.

Fig. 11: Input and output in a gate region. The boundary is divided into 6 segments, which are in different color. Blue arrows represent possible outputs. Orange arrows represent possible inputs.

For a discrete-time Markov process, the continuous states needed to be discretized into a finite number of input and output states. Here, the state of the body when crossing the boundary is , where {EN, RD, LD, RT, LT, MT} is the boundary index, is the position on the corresponding the boundary. is measured as the distance to several benchmark points (EN: m, LT/LD: m, RT/RD: m, MT: m). is the velocity when crossing boundaries. To set up the state space for the Markov chain, we discretized the input states for boundaries EN, RD and LD (Table III). Inputs from RT or LT rarely occur (fewer than 4% of all trials) and they are not involved with the body-beam interaction. Thus, for simplicity, we did not further discrete RT and LT like the other boundaries to improve computational efficiency. In total, there are 87 input states, , and 88 output states, , with being the state of being trapped in a gate.

[-25,-15], [-15,-5],
[-5,5], [5,15],
[10,20],[20,30] [10,20],[20,30]
[-25,-15], [-15,-5],
[-5,5], [5,15],
(m/s) [10,20]
Groups 25 30 30
TABLE III: Discretization of Input States

For Markov Chain analysis, we need to obtain the transition matrix , where is the probability from input to the output during a gate traversal (Eqn. 16). For our discretization, is an matrix. The matrix is obtained by running 100 trials per group with random force in the single gate simulation.


where is the probability of output state after the traversal, is probability of input state before the traversal.

With the transition matrix, the Markov chain Monte Carlo (MCMC) method [8]

can be used to predict the consecutive transitions between gate units. Monte Carlo is a statistical simulation method relying on repeated random sampling from a probability distribution. MCMC is the combination of Markov chain and Monte Carlo. In each gate, input is known, either from the initialization or from the former gate. The probability of outputs after traversing this gate is given by multiplying the Markov chain transition matrix. We randomly pick a certain output state based on the probability distribution

and use it as the input of the next gate. By repeating the Markov chain and Monte Carlo in turn, we can obtain a path across the terrain. The simulation ends after a maximal number of iterations. The pseudo code:

while Body within terrain do
     Markov Chain
     Randomize one output Sample from
     if = then Break
     end if
     Transition to the next gate, new input
     Reset zeros(1, 87)
end while
Algorithm Markov Chain Monte Carlo

Unlike our dynamic simulation, the Markov chain Monte Carlo algorithm cannot precisely predict the trajectory across the terrain. But we can compare the prediction for the finial lattices where the trials finally stop. By repeating MCMC trials, we can obtain a predicted distribution of the final lattice. To compare two methods, we have the same parameters of gates in dynamic simulation as the ones in MCMC. Ideally, the MCMC simulation can predict the path on an infinitely large terrain. However, if the body keeps accelerating after traversing a series of gates, the locomotor velocity will keep increasing towards infinity, and it would be impractical to discretize velocity for the Markov chain analysis. To address this issue, we increased the beam damping ( Ns/rad) and added a viscous force ( Ns/m) between body and the horizontal plane so that the velocity can remain in the tested range (Table. III). The body started at the origin of the entire terrain ( = 0 m, = 0 m) with = 0 m/s and = 15 m/s. For our gate lattice, each MCMC trial is allowed to have a maximum of 13 steps to reach the furthest gate. But if the body exits from the terrain or becomes trapped inside a gate (i.e., reaching ), the trial will stop at that gate.

To test how well the MCMC simulation predicted traversal outcome, we ran 100 trials using both dynamic and MCMC simulation and compared the distribution of the final location of the body within the obstacle field. The patterns are strikingly similar (Fig. 12). The correlation coefficient, , between the two simulations was 0.914. In addition, the root mean squared error, , was only 1.018, i.e., the error is around 1 out of 100 trials. These results demonstrated that the MCMC simulation can predict our dynamic simulation traversal result very well. Considering that the MCMC method can finish the 100-trial prediction within 10 seconds whereas the dynamics simulation took around 100 hours, there is a huge advantage in saving computational time of the MCMC method.

Fig. 12: Distribution of the final location of the body within the obstacle field. (A) The predicted distribution given by MCMC simulation (left) and dynamic simulation (right). (B) Data in A, rearranged into a more direct pair-wise comparison.

To evaluate the generality of this result, we tested three other sets of randomly chosen initial states and found excellent match throughout (Table IV).

(m) (m/s) (m/s)
0 0 15 0.914 1.018
23 11 16 0.9097 1.133
4 5 10 0.8929 1.252
17 13 12 0.9026 0.9813
TABLE IV: Prediction Evaluation Using Different Initial States.

Iv Discussion

In this study, we developed a stochastic dynamic model to simulate a self-propelled, minimalistic locomotor body with random forces traversing beam obstacles in two dimensions on a horizontal plane. We found that a larger forward propulsive force and a sufficient but not excessive lateral random force facilitate traversal, in accord with the kinetic energy needed to overcome potential energy barriers. Asymmetry because of different horizontal beams stiffness can cause asymmetric distribution in trajectories, in accord with the higher probability to traverse on the lower barrier side. The model of interacting with a basic unit of obstacle (two beams) was further extended to simulate the traversal of a multi-obstacle field. Finally, we applied the Markov chain Monte Carlo method to predict how the body statistically transitions between adjacent gates, which predicted stochastic paths on the larger terrain well.

Our results demonstrated that the potential energy landscape is useful for understanding traversal dynamics without solving them exactly. This would be particularly useful for modeling systems with complex interaction, which is common in animal and robot locomotion in complex terrain. In addition, we can redesign the landscape barriers and basins by adjusting the model parameters related to physical interaction, which can be applied in robotics challenges such as motion control and path planning. Furthermore, stochasticity, which is common in locomotion and usually considered a nuisance, can be beneficial for robots when they are trapped by obstacles and need to overcome potential energy barriers. Finally, the exploration of multi-obstacle field traversal supported the plausibility to decompose natural terrain as a combination of multiple obstacles [16]. Our results suggested that, by statistically understanding and predicting the interaction with each type of obstacle unit, larger scale traversal processes may be rapidly predicted.


The authors would like to thank Noah Cowan, Ratan Othayoth, and Yaqing Wang for discussion and three anonymous reviewers for suggestions.


  • [1] A. D. Ames, A. Abate, and S. Sastry (2005) Sufficient conditions for the existence of Zeno behavior. In Proceedings of the 44th IEEE Conference on Decision and Control, pp. 696–701. Cited by: §II-B.
  • [2] C. M. Bergman, J. A. Schaefer, and S. Luttich (2000) Caribou movement as a correlated random walk. Oecologia 123 (3), pp. 364–374. Cited by: §I.
  • [3] J. Borenstein, Y. Koren, et al. (1991)

    The vector field histogram-fast obstacle avoidance for mobile robots

    IEEE Transactions on Robotics and Automation 7 (3), pp. 278–288. Cited by: §I.
  • [4] G. N. DeSouza and A. C. Kak (2002) Vision for mobile robot navigation: a survey. IEEE Transactions on Pattern Analysis and Machine Intelligence 24 (2), pp. 237–267. Cited by: §I.
  • [5] A. A. Faisal, L. P. Selen, and D. M. Wolpert (2008) Noise in the nervous system. Nature Reviews Neuroscience 9 (4), pp. 292–303. Cited by: §I.
  • [6] S. W. Gart and C. Li (2018) Body-terrain interaction affects large bump traversal of insects and legged robots. Bioinspiration & biomimetics 13 (2), pp. 026005. Cited by: §I.
  • [7] S. W. Gart, C. Yan, R. Othayoth, Z. Ren, and C. Li (2018) Dynamic traversal of large gaps by insects and legged robots reveals a template. Bioinspiration & biomimetics 13 (2), pp. 026006. Cited by: §I.
  • [8] W. R. Gilks, S. Richardson, and D. Spiegelhalter (1995) Markov chain monte carlo in practice. CRC press. Cited by: §III-C.
  • [9] Y. Han, R. Othayoth, Y. Wang, C. Hsu, R. de la Tijera Obert, E. Francois, and C. Li (2021) Shape-induced obstacle attraction and repulsion during dynamic locomotion. The International Journal of Robotics Research 40 (6-7), pp. 939–955. Cited by: §I, §I.
  • [10] C. Harley, B. English, and R. Ritzmann (2009) Characterization of obstacle negotiation behaviors in the cockroach, Blaberus discoidalis. Journal of Experimental Biology 212 (10), pp. 1463–1476. Cited by: §I.
  • [11] C. Li, A. O. Pullin, D. W. Haldane, H. K. Lam, R. S. Fearing, and R. J. Full (2015) Terradynamically streamlined shapes in animals and robots enhance traversability through densely cluttered terrain. Bioinspiration & Biomimetics 10 (4), pp. 046003. Cited by: §I.
  • [12] G. Li, W. Li, J. Zhang, and H. Zhang (2015) Analysis and design of asymmetric oscillation for caterpillar-like locomotion. Journal of Bionic Engineering 12 (2), pp. 190–203. Cited by: §I.
  • [13] T. Y. Moore, K. L. Cooper, A. A. Biewener, and R. Vasudevan (2017) Unpredictability of escape trajectory explains predator evasion ability and microhabitat preference of desert rodents. Nature Communications 8 (1), pp. 1–9. Cited by: §I, §II-A.
  • [14] J. R. Norris (1998) Markov chains. Cambridge university press. Cited by: §III-C.
  • [15] R. Othayoth, G. Thoms, and C. Li (2020) An energy landscape approach to locomotor transitions in complex 3d terrain. Proceedings of the National Academy of Sciences 117 (26), pp. 14987–14995. Cited by: §I, §I, §III-A.
  • [16] R. Othayoth, Q. Xuan, Y. Wang, and C. Li (2021) Locomotor transitions in the potential energy landscape-dominated regime. Proceedings of the Royal Society B: Biological Sciences 288 (1949). Cited by: §I, §I, §I, §II-C, §IV.
  • [17] E. Rimon (1990) Exact robot navigation using artificial potential functions. Ph.D. Thesis, Yale University. Cited by: §I.
  • [18] J. Schmitt and P. Holmes (2000) Mechanical models for insect locomotion: dynamics and stability in the horizontal plane I. Theory. Biological Cybernetics 83 (6), pp. 501–515. Cited by: §II-A.
  • [19] B. P. Selgrade, M. Thajchayapong, G. E. Lee, M. E. Toney, and Y. Chang (2017) Changes in mechanical work during neural adaptation to asymmetric locomotion. Journal of Experimental Biology 220 (16), pp. 2993–3000. Cited by: §I.
  • [20] Q. Xuan and C. Li (2020) Coordinated appendages accumulate more energy to self-right on the ground. IEEE Robotics and Automation Letters 5 (4), pp. 6137–6144. Cited by: §I.
  • [21] Q. Xuan and C. Li (2020) Randomness in appendage coordination facilitates strenuous ground self-righting. Bioinspiration & Biomimetics 15 (6), pp. 065004. Cited by: §I.