Environmental force sensing helps robots traverse cluttered large obstacles

Robots can traverse sparse obstacles by sensing environmental geometry and avoiding contact with obstacles. However, for search and rescue in rubble, environmental monitoring through dense vegetation, and planetary exploration over Martian and lunar rocks, robots must traverse cluttered obstacles as large as themselves by physically interacting with them. Previous work discovered that the forest floor-dwelling discoid cockroach and a sensor-less minimalistic robot can traverse cluttered grass-like beam obstacles of various stiffness by transitioning across different locomotor modes. Yet the animal was better at traversal than the sensor-less robot, likely by sensing forces during obstacle interaction to control its locomotor transitions. Inspired by this, here we demonstrated in simulation that environmental force sensing helps robots traverse cluttered large obstacles. First, we developed a multi-body dynamics simulation and a physics model of the minimalistic robot interacting with beams to estimate beam stiffness from the sensed contact forces. Then, we developed a force feedback strategy for the robot to use the sensed beam stiffness to choose the locomotor mode with a lower mechanical energy cost. With feedforward pushing, the robot was stuck in front of stiff beams if it has a limited force capacity; without force limit, it traversed but suffered a high energy cost. Using obstacle avoidance, the robot traversed beams by avoiding beam contact regardless of beam stiffness, resulting in a high energy cost for flimsy beams. With force feedback, the robot determined beam stiffness, then traversed flimsy beams by pushing them over and stiff beams by rolling through the gap between them with a low energy cost. Stiffness estimation based on force sensing was accurate across varied body oscillation amplitude and frequency and position sensing uncertainty. Mechanical energy cost of traversal increased with sensorimotor delay. Future work should demonstrate cluttered large obstacle traversal using force feedback in a physical robot.


Introduction
Mobile robots are becoming increasingly prevalent in society.Many robots already excel at navigating on flat surface and avoiding sparse obstacles by sensing environmental geometry, like robot vacuum cleaners [1] and self-driving vehicles [2].Sensors like cameras, radars, and LiDAR are used to create a geometric map of the surroundings, then a collision-free path is planned to avoid physical contact with obstacles [3][4][5][6][7].However, for search and rescue through earthquake rubble [8,9], environmental monitoring through dense vegetation [10], and extraterrestrial exploration through Martian and lunar rocks [11], robots often need to traverse complex 3D terrain with densely cluttered obstacles as large as themselves by physically interact with the obstacles.In these situations, geometry-sensing based obstacle avoidance faces difficulties.Firstly, a collision-free path might not always be available.Secondly, even if such a path exists, it could be energetically costly or need an extended execution time, as the robot may need to take long detours to circumvent obstacles.Additionally, if the robot fails to avoid obstacles promptly, unexpected collisions and contact could cause it to flip over [12,13].
Understanding how to utilize and control physical interaction with the environment (environmental affordance [14,15]) is fundamental to locomotion.For example, understanding of terramechanics [16,17] has helped wheeled vehicles off-road deformable ground.Understanding leg-ground interaction on rigid ground [18][19][20] has provided the foundation for legged robots to walk and run on flat surfaces with small unevenness.Recent research into leg interaction with granular media has facilitated the robotic design and control to improve legged and even (surprisingly) wheeled locomotion on deformable ground [21][22][23][24][25][26].
Biological organisms provide rich sources of inspiration for understanding and controlling the physical interaction of robots with obstacles to improve their traversal in complex 3D terrains with cluttered large obstacles.An excellent model organism is the forest floor-dwelling discoid cockroach recognized for its exceptional ability to traverse vegetation, foliage, and rocky crevices [27,28].Recent studies have shed considerable light on how these creatures use effective body and leg interaction with obstacles to traverse cluttered large obstacles via integration of biological experimentation [27,29,30], robophysical modeling [29,[31][32][33][34][35], and theoretical and computational modeling [36].
Prior research discovered that the discoid cockroach can traverse grass-like beam obstacles using and transitioning across different modes of locomotion [27].For example, the animal simply pitches up the body to push over flimsy beams (figure 1(A), blue), whereas it often transitions (figure 1(A), yellow) to rolling (figure 1(A), red) into a gap to maneuver through stiff beams [35].Its legs continuously pushed against the ground, inducing oscillatory selfpropulsion during traversal [35].
To study the principles of locomotor mode transitions, a sensor-less, feedforward, minimalistic robophysical model was developed (figure 1(B)) [35].The robot, featuring an ellipsoid-like body without legs, was propelled forward at a consistent speed.Its body oscillated to emulate the oscillatory selfpropulsion of the animal.This robot was free to pitch and roll in response to interaction with two beams.The grass-like beams were modeled as two plates with torsion springs at their base on the ground.Utilizing this robot, researchers discovered the mechanism of locomotor mode transitions under passive interaction with obstacles [35].A potential energy landscape model showed that the pitch and roll modes arose as the system was attracted to distinct pitch and roll basins on a potential energy landscape, respectively [35].Systematic experiments using this robot demonstrated that locomotor transition from pitch to roll mode occurred when kinetic energy fluctuation from body oscillations exceeded the potential energy barrier between the pitch and roll basins [35].
However, although both the animal and robot can traverse cluttered grass-like beam obstacles of a range of stiffness, the animal was still better at traversing.The animal often made the transition even when the kinetic energy fluctuation from body oscillations alone was insufficient for overcoming the potential energy barrier [35], possibly by utilizing force sensing during environmental interaction to actively control its locomotor mode transitions.This was partly substantiated by a subsequent study which found that the animal actively adjusted their body and appendages (e.g.head, abdomen, and legs) when transitioning from the pitch to the roll mode [30].
Here, we further studied how force sensing during body-obstacle interaction could help active locomotor mode transitions.We chose to use the minimalistic, legless robot design because it provides an ideal, simplified system to focus on sensing and controlling body-obstacle interaction, which is important to locomotor transitions through cluttered large obstacles [35].The potential energy landscape model in the previous study did not fully model the dynamics of the system [35], limiting our ability to explore the impact of force sensing.To address this, we developed a multi-body dynamics simulation of the minimalistic robot with force sensing and feedback (section 2.1).While environmental force sensing capabilities have been developed for some bio-inspired legged robots [37][38][39][40][41][42], they only detect contact force between the robot feet and the ground.Creating a refined robot capable of sensing body-obstacle contact forces and torques requires integrating distributed force sensors that cover the body surface and developing feedback control algorithms to control body rotations.This is a major undertaking and is beyond the scope of this study.
We developed a force feedback strategy to use the robot's body force sensing to estimate the physical properties of the beams and alter the traversal method for different beams.In addition, we developed a physics model of the system that integrated a contact force model and the corresponding equations of motion (section 2.2).The physics model allowed us to estimate the beam stiffness from the contact forces.We performed simulation studies to compare strategies with and without force sensing (section 2.3).Our force feedback strategy first estimated beam stiffness based on force sensing and a contact force model (section 2.4).We then used the potential energy landscape model to classify beams as either stiff or flimsy to guide the robot to choose a locomotor mode with a lower energy cost (section 2.5).Using this information, we performed motion planning to generate robot trajectory (section 2.6).Our simulation also allowed us to vary body oscillations and randomness during traversal to examine their influence on Discoid cockroach and robophysical model traverse grass-like beam obstacles using physical interaction.(A) For flimsy beams, the discoid cockroach predominantly uses a pitch mode (blue), while for stiff beams, it predominantly transitions (yellow) to a roll mode (red) [35].(B) A minimalistic robophysical model of cockroach beam traversal [35].Adapted with permission from [35].our strategy.Oscillation is a naturally occurring phenomenon during traversal.The discoid cockroach displays substantial body oscillations while traversing due to their legs periodically pushing against the ground [35].Oscillations lead to intermittent contact instead of steady contact, which can make position sensing less accurate and presents increased challenges in modeling contact forces and estimating stiffness.
To evaluate the usefulness of the force feedback strategy, we conducted a series of simulation studies.First, we analyzed simulations of a robot traversing both stiff and flimsy beams without force sensing and studied its three different strategies (section 3.1).We studied feedforward pushing strategy (section 3.1.1)and obstacle avoidance strategy (section 3.1.3)to provide a baseline.Considering that in real-world scenarios both animals and robots have finite force/torque capacities, we next studied feedforward pushing strategy with force/torque limits (section 3.1.2).Second, we tested the robustness of the force feedback strategy against body oscillations with randomness and inaccuracies in position sensing (section 3.2) to compare with the baseline.Third, to access the advantages of the force sensing and active feedback control, we conducted simulations with the force feedback control strategy and compared results with those from previous non-force sensing strategies (section 3.3).Furthermore, as the sensorimotor delay is inevitable in both animals and robots, we used our simulation to study the influence of sensorimotor delay on traversal performance (section 3.4).Finally, we discussed the limitations of our method, its broader applications, and future work (section 4).

Methods
In this section, we first describe how we developed a physics-based simulation of a minimalistic robot interacting with beams, incorporating body force sensing (section 2.1).Then, we define the physics model that integrated a contact force model along with the corresponding equations of motion (section 2.2).Next, we describe how the force feedback strategy used force sensing and potential energy landscape to guide the robot to select a locomotor mode that minimized energy consumption (sections 2.4, 2.5 and 2.6).We formulated the energy cost during traversal and conducted simulations to compare strategies with and without force sensing (sections 2.3, 2.7 and 2.8).Lastly, we experimented with varying the sensorimotor delay in force feedback control (section 2.9).

Simulation of robot traversing beams with body force sensing
We developed a physics-based simulation as a platform to explore the benefits of sensing contact forces between the robot body and obstacles and using this information to enhance its ability to traverse beams.The simulation robot (figure 2(A)) was built in Chrono, an open-source multi-body dynamics engine [43,44], following the design in previous study [35].It was modeled as a rigid ellipsoidal body similar to the discoid cockroach's body, with principal axes lengths 2a = 0.22 m, 2b = 0.16 m, and 2 c = 0.06 m.An additional weight of 0.5 kg made it bottom-heavy, with the center of mass at h c = 0.01 m below the geometric center (figure 3(A)).The ellipsoidal body's total mass M = 1 kg and its moments of inertia along the three principal axes were (I x , I y , I z ) = (3.9,7.5, 10.9) × 10 −3 kg m 2 .The geometric center was H o = 0.105 m from the ground before moving.
The ellipsoidal body had three degrees of freedom in translation and two degrees of freedom in rotation (figure 2(A)).Although the robot was equipped with three linear actuators for 3D translation, it was primarily propelled forward by a linear actuator along the X-axis.The lateral and vertical actuators only generated oscillations that allow for slight movements (∼2 mm) along the Y-and Z-axes.In accordance with the previous design [35], two rotation axes, which were the roll axis (figure 2, red dashed line) and the pitch axis (figure 2, blue dashed line).Two motors controlled the ellipsoidal body's rotation to emulate the animal's active pitch and roll adjustments.The center of the rotation overlapped with the geometric center.Thus, body orientation can be represented by roll angle α and pitch angle β.The body was in static equilibrium when both its roll angle α and pitch angle β were 0 • .Following previous work [35], the grass-like beams were designed as rigid, vertical plates with torsion springs at their base, and they could deflect about Y-axis.The torque function of the spring was where k was stiffness, θ was the relative beam deflection angle, c d was the damping coefficient (c d = 0.01 Nm s/rad), and ω was the beam's deflection angular velocity.Note that here we added an additional damping term to represent the energy dissipation occurring during the beam's motion.The mass of each beam was m = 1 g.The two torsion springs were symmetric about the X-Z plane.The two beams were vertical when there was no deflection.The gap width between the two beams was d = 0.138 m, which was narrower than the robot body width.Each beam's height, width, and thickness were L = 0.155 m, w = 0.04 m, and h = 0.01 m (figure 3(B)).
We used the discrete-element method with penalty in the Chrono framework, which characterizes contact through a viscoelastic force model [44,45] . (1) Here, the vector sum of δ n (normal) and δ t (tangential) was the displacement vector at the contact point between two bodies, which symbolizes their overlap in the absence of deformation.The effective mass m and the effective radius of curvature R are defined as m = m1m2 m1+m2 and R = R1R2 R1+R2 .k n , k t , γ n , γ t are the normal and tangential stiffness and damping parameters, respectively.They are all dependent on m, R, and the material properties of the interacting bodies, including Young's modulus (E), Poisson's ratio (ν), and the coefficient of restitution (CoR).In our particular design, contact occurred exclusively between the robot's body and the beams.We set their Young's modulus = 10 6 Pa, Poison's ratio = 0.1, coefficient of restitution = 0, and coefficient of friction = 0. Within the Chrono simulation, we could determine the position and orientation of the body and beams (environmental geometry sensing), as well as their contact forces (environmental force sensing).Thus, in simulation, we could control the motors based on these sensory inputs (feedback control).
To validate the simulation, we replicated the experiments from previous research conducted with a real robot [35].Consistent with the previous findings, we verified that the oscillatory motion aided the robot in transitioning from pitching to rolling to traverse beams via passive interaction with the beams (figures 2(B) and (C)).

Physics model of body-beams interaction
We built a physics model (figure 3) based on the simulation robot design to estimate the beams' stiffness using environmental force sensing.The physics model consisted of the contact force model and equations of motion.
Besides the input force F x and torques (τ 1 , τ 2 ) from actuators acting on the simulation robot (figure 3(A)), there were three external forcesgravitational force ⃗ G and contact forces ( ⃗ F 1 , ⃗ F 2 ) with two beams (figure 3(C)).Given these, we could calculate the total force ⃗ F and torque ⃗ τ .Thus, applying Newton-Euler equations, the equations of motion of the robot were: where M was the mass, I CoM was the moment of inertia about the center of mass, ⃗ α and ⃗ β were acceleration and angular acceleration about the center of mass, and ⃗ ω was the angular velocity.In addition, we assumed that actuators would not recycle energy.For example, when the body rolls back to a lower energy state, the decrease in mechanical energy is not returned.Thus, energy is dissipated by motor doing negative work, leading to a mechanical energy cost.
Because the beams were much lighter than the robot, the inertial force of the beam was negligible compared to the contact forces.Therefore, we assumed that each beam was quasi-static and satisfied torque equilibrium during contact where k i and c di were the stiffness and damping coefficients of the ith beam, θ i and ω i were its beam deflection angle and angular velocity,⃗ r bi was the distance from the revolute axis of the ith beam to the point of application of ⃗ F i .
Because the robot's body surface was smooth, the coefficient of friction was set to zero in the simulation.Thus, the contact force was normal to the tangent plane of the contact point on the ellipsoidal body's surface.
The normal direction of the ellipsoidal body's tangent plane at a point (x, y, z) in the body frame was: Given the current position and orientation of the robot and equations ( 3) and ( 4), the theoretical contact forces ⃗ F 1 and ⃗ F 2 in the world frame was a function of stiffness k 1 and k 2 .

Simulation studies to compare between without and with force sensing
We conducted simulations to investigate the traversal of beams with varying stiffness under different strategies.We first studied traversing both flimsy and stiff beams (defined in section 2.5) without force sensing, using three different strategies: feedforward pushing without force/torque limits (section 3.1.1),feedforward pushing with force/torque limits (section 3.1.2),and avoiding obstacles (section 3.1.3).Then, we studied traversal of both stiff and flimsy beams using our force feedback strategy (section 3.3).
For the feedforward pushing strategy without force/torque limits, the robot had an arbitrarily large propulsive force to maintain a forward speed of v X = 0.05 m s −1 , with no control exerted over pitch and roll directions.
For the feedforward pushing strategy with force/torque limits, we set 1 N m to examine more realistic traversal with force and limits.When any force or torque component of ⃗ u calculated from equation ( 7) exceeded its limit, this component was set to be the maximum or minimum values within the limit.
For the obstacle avoidance strategy, we used environmental geometry to determine the minimal angle required for the robot to roll to traverse the beams without contacting them.The robot maintained a forward speed of v X = 0.05 m s −1 utilizing roll control to rotate its body and avoid contact.Although the above force/torque limits were also applied in this scenario, the absence of contact resulted in no resistance, so these limits were not reached.
For the force feedback strategy, we first estimated beam stiffness based on the body contact force (sensing time T s = 100 ms) sensed after the initial contact.Then, we constructed a potential energy landscape and performed motion planning.Finally, we controlled the robot to track planned trajectories to traverse flimsy and stiff beams.Detailed discussions on the motion planning and control were provided in sections 2.6 and 2.7.To be fair in comparison, the same force/torque limits were also applied in this scenario.

Estimating beam stiffness from force sensing
The vector sum of contact forces between the robot and two beams were sensed with a sampling rate of 40 Hz, which were used to calculate the unknown beam stiffness k 1 and k 2 in equation ( 3) by minimizing the error between the theoretical contact force ⃗ F = ⃗ F 1 + ⃗ F 2 and the measured force data ⃗ F sensed (using the MATLAB 'fminsearch' function): where N was the total number of the sampled force data points at t i for i = 1 to N. To account for the fact that this method could only determine a local minimum, this process was repeated 100 times using various randomly chosen initial guesses for stiffness, increasing the likelihood of finding a global minimum.Initial guess values of k 1, 2 were chosen randomly from a uniform distribution within the range of [0, 5] N m rad −1 , which was the range of interest for our study.
To quantify the estimation accuracy, we calculated the relative error of stiffness, which was the ratio between the estimated error and its true value where k i is the estimate value (i = 1, 2), k i is the true value of stiffness (k 1 = k 2 = 0.5 N m rad −1 ).
To study the robustness of the estimation accuracy, we studied the effect of random body oscillations, inaccurate position sensing, and sensing time.

Varying vertical and lateral body oscillations with randomness
To investigate the impact of naturally occurring random body oscillations [35] on our beam stiffness estimation method, we simulated the robot's movement with lateral and vertical oscillations, generated by linear actuators (figure 2(A)).The lateral oscillation followed the triangle wave function, with the frequency f and amplitude A. The vertical oscillation followed a sum of 30 sine functions to induce randomness.The sum of sine functions is a strategy employed in various research contexts to introduce unpredictable external stimuli [46][47][48].The amplitude of each sine function was randomly selected from −0.5 mm to 0.5 mm, and the phase was randomly chosen from −π to π.The frequencies of 30 sine functions were f 50 , 2f 50 , 3f 50 , …, 30f 50 .We conducted the simulations with varying oscillation amplitude (A = 1, 2, 3 mm) and frequency (f = 2, 4, 6 Hz), while assuming that the position sensing remained accurate.We ran five trials for each set of oscillation amplitude and frequency.

Varying inaccuracy of position sensing
To study the impact of uncertainty in the sensed position on the accuracy of beam stiffness estimation, we used an inaccurate sensed Y position, by assuming Y = 0 mm during traversal regardless of lateral oscillation.Thus, the larger the lateral oscillation amplitude, the less accurate the sensed Y position was.We ran simulations with an oscillation frequency of 6 Hz and varied the oscillation amplitude = 1, 2, and 3 mm.For each oscillation amplitude, we ran five trials.

Varying sensing time
Given our constant forward speed, it only took the robot 4.4 s to traverse a distance equal to its body length.Thus, the robot needs to react fast, i.e., the time cost of our method is important.
An inevitable time cost is sensing time, which is required to obtain contact force data to estimate stiffness.Therefore, we studied how beam stiffness estimation accuracy depends on the sensing time, to inform how long the robot needs to sense to obtain a reasonable estimate.We conducted tests using sensing time T s = 25, 50, 100, and 200 ms, with a data sampling rate of 40 Hz.These tests were collected at an oscillation amplitude of 1 mm and a frequency of 2 Hz, with five trials executed for each individual sensing time.

Potential energy barrier calculation to differentiate between stiff and flimsy beams
Cockroaches use different strategies when traversing flimsy and stiff beams [35].A critical stiffness helps categorize whether beams are stiff and should be traversed via the role mode (rolling through the gap between beams), or whether they are flimsy and should be traversed via the pitch mode (push over the beams).
To find the critical stiffness, we calculated the potential energy barrier corresponding to the pitch and roll modes.This allowed us to estimate which mode costs less mechanical energy for traversing stiff or flimsy beams.The potential energy of the system consisted of the gravitational potential energy of the robot and the elastic potential energy of the beams: where Z CoM (α, β) was the height of the center of mass from the ground, θ 1 and θ 2 were beam deflection angles, which were functions of α, β, and X.We varied X from −0.12 m to 0.12 m with an increment of ∆X = 0.002 m and α and β from −90 • to 90 • with an increment of ∆α or ∆β = 2 • to obtain a 3D potential energy landscape over the entire (α, β, X) workspace during traversal (figure 4(A)).
For simplicity, we started with the case where two beams had the same stiffness (k 1 = k 2 = k 0 ).
For the pitch mode, we assumed that the potential energy reaches the maximum when the robot's bottom reaches the top of the beams where beams were deflected most (figure 3(D), left).Thus, its potential energy barrier was PE pitch (θ) = 2 × 1 2 k 0 θ 2 , which was the elastic potential energy of two beams.Because cos θ 0 = H L , we had PE pitch = k 0 (cos −1 H L ) 2 .For the roll mode, the potential energy reached the maximum when the robot rolled until there was no contact in the entire traversal (figure 3(D), right).Thus, its potential energy barrier was PE roll = Mg • ∆h, which was the increase of gravitational potential energy of the body.Using the principles of similar triangles, we obtained ∆h = h c 1 − d 2b , which was the lifted height of center of mass of the body (figure 3

(D), right).
When PE pitch − PE roll = 0, the threshold beam . Because PE roll did not depend on beam stiffness whereas PE pitch was proportional to it, the barrier difference PE pitch − PE roll monotonically increased with k 0 .Above this threshold, the roll mode was advantageous due to its lower energy barrier, whereas below this threshold, the pitch mode had a lower barrier.Substituting the parameters' values given in the simulation robot design (section 2.1), we obtained k 0 = 0.146 N m rad −1 .
For simplicity, we selected a single representative stiffness value for each category: k low = 0.01 N m rad −1 and k high = 0.2 N m rad −1 as the stiffness of flimsy and stiff beams, respectively.

Force feedback strategy motion planning
For the force feedback strategy, we need to obtain the desired trajectory for traversing beams in X, roll, and pitch space according to beams stiffness.To determine this trajectory, a prediction of the mechanical energy cost associated with a given path was necessary.A weighted directed graph (figure 4(B)), composed of vertices and directed edges, was constructed based on the energy landscape (section 2.5).Each vertex was connected through directed edges to its six closest neighbors, (X ± ∆X, α, β), (X, α ± ∆α, β), and (X, α, β ± ∆β), where the increments were defined in section 2.5.The directed edge had a value representing the mechanical energy cost.The total mechanical energy cost of a path in the network equaled the sum of the values attributed to its constituent directed edges.In a straightforward prediction and approximation, we defined the value of a directed edge leading to a state with greater potential energy to be positive, quantified by the increase in this potential energy.The value of a directed edge heading towards a state with lower potential energy was defined as zero, denoting no energy cost and recovery for such transitions (figure 4(B)).
For each trajectory, we defined the initial state as the state when the robot started the active control (t s ), following a sensorimotor delay ∆t from its first contact with the beams (t c ).The sensorimotor delay was ∆t = t s − t c .The target state was defined as the state after traversal at X = 0.1 m with the horizontal body pose (α = β = 0 • ).Given the directed weighed graph (figure 4(B)), the A-star algorithm [49] was used to search for the path of the minimal cost between the initial and target states (figure 4(C)), with the heuristic function being the potential energy difference of two states.The A-star algorithm operates by evaluating the cost to travel from one node to another, coupled with a heuristic that estimates the cost to reach the goal from a given node.We took this path as the planned path.

Force feedback strategy motion control
For the force feedback strategy, motion control was required to track the planned trajectory.We controlled the roll and pitch angle while maintaining a forward speed of v X = 0.05 m s −1 .This constant forward speed setting simplified the control process and was retained from the prior design [35].Based on the position error, we used a model-based feedback control to obtain the control input at every 0.002 s: where ⃗ q = (X, α, β) was the current state and − → q d = (X d , α d , β d ) was the desired state.K p was positive feedback gain (we chose 0.1).⃗ u = (τ 1 , τ 2 , F X ) was the control input (figure 3(B)).τ 1 , τ 2 are pitch and roll torques.F X is the force along X direction.F ext (⃗ q) was the external force calculated from the contact force model (section 2.2), including contact forces with beams and the gravitational force.With the control input ⃗ u, the current state was controlled to approach the desired state at every time step.

Motor work calculation to compare across strategies
To compare the mechanical energy cost of different strategies, we calculated the work done by actuators in the entire traversal using the below function: where m was the number of control steps in the entire traversal time, dt = 0.002 s was the control time step (section 2.7).(F X v X ) was the power of the force F X .

Varying sensorimotor delay for force feedback strategy
Sensorimotor delay is inevitable in both animals and robots.To study the impact of sensorimotor delay on traversal performance of our force feedback strategy, we conducted simulations with three different sensorimotor delays ∆t = 320, 480, and 640 ms and evaluated its impact on mechanical energy cost from actuators (equation ( 8)).Sensorimotor delay was defined as the time interval between the instant of the first contact with beams and the start of actuators, incorporating the previously mentioned sensing time T s .
For simplicity, and to better show the impact of sensorimotor delay, the initial body pitch angle started with zero angle and the actuators did not generate body oscillations.Because the simulation was deterministic without randomness, we ran one trial for each sensorimotor delay.

Results
We found that each of the three strategies without force sensing offers distinct advantages and drawbacks when traversing flimsy and stiff beams (section 3.1).Our force feedback strategy had the combined strengths of these methods, facilitating enhanced beam traversal of both flimsy and stiff beams (section 3.3).Beam stiff estimation was robust against random oscillations (section 3.2).Additionally, the mechanical energy cost during traversal increased with sensorimotor delay for the force feedback strategy (section 3.4).

Traversing flimsy and stiff beams without force sensing 3.1.1. Feedforward pushing without force/torque limits
When encountering flimsy beams, the robot pushed over beams with a small body pitch angle (<4 • ) (figure 5(A), bottom; movie 1) and small contact force (<0.13 N).Its energy cost over the entire traversal was 14.8 mJ.When encountering stiff beams, the robot was pushed by the beams to oscillate up and down over a wide range of angles ([−37 • , 38 • ]) (figure 5(A), top; movie 1), and the maximal resistive force from the beams reached 1.5 N. If the robot were freely moving on the ground and experienced such large pitch oscillations, it could easily lose stability and flip over.In addition, the robot's mechanical energy cost over the entire traversal was 185.2 mJ, much larger than during flimsy beam traversal.

Feedforward pushing with force/torque limits
When encountering flimsy beams, the robot behaved the same as it did with unlimited force because the resistive force was smaller than the force limit (figure 5(B), bottom; movie 1).However, when encountering stiff beams, the robot was stuck in front of the beams (figure 5(B), top; movie 1).The body pitched up substantially in this process, and it then oscillated up and down in pitch from intermittent contact with beams.

Avoiding obstacles
Because beams with different stiffness had the same environmental geometry, the robot moved in the exact same way when traversing flimsy and stiff beams (figure 5(C); movie 1).Its energy cost was 35.8 mJ.The obstacle avoidance strategy did save much energy compared to the feedforward pushing strategy (185.2 mJ) when the obstacle was stiff.However, the obstacle avoidance strategy was not the best option for the flimsy beams.Its energy cost was higher than the feedforward pushing strategy (14.8 mJ), and the large body rolling required to avoid the beams could make it challenging to maintain stability.

Estimating beam properties in simulation
We estimated the stiffness of the beams using the contact physics model and studied the robustness of the estimation accuracy against random body oscillations, inaccuracy position sensing, and the varying sensing time.

Estimation of beam stiffness under random body oscillations
Our estimation results for two beams' stiffness ( k 1 , k 2 ) across varying oscillation frequencies indicated that the estimation accuracy was not notably affected by the oscillation frequencies (table 1) or amplitudes (table 2).Although the oscillation caused more collisions and intermittent body-beam contact, our approach still obtained an accurate estimation of stiffness (relative error e k < 5% in 90% trials), which was a prerequisite for motion planning and control.

Estimation of beam stiffness with position sensing error
In the above estimation, although oscillations and randomness were introduced, the robot still had accurate measurements of its position relative to the obstacles.However, in practice, the robot's position sensing may not be accurate because of its intense oscillations.
Through calculating the mean of relative errors ( e k1 +e k2 2 ) of beams' stiffness among trials, we discovered that, generally, an increase in oscillation amplitude led to an increase in the relative error of estimation (figure 6) as well as an increase in the variation of estimation among different trials.Despite these, the estimation accuracy was still acceptable (relative error e k < 15%), which was sufficient for differentiating stiff and flimsy beams.Thus, our estimation was robust against errors in position sensing.).This sum represents the overall error.However, for a sufficiently long sensing time (⩾50 ms), further increasing sensing time (adding more sensed data) did not improve estimation accuracy substantially, because the randomness and modeling error became dominant (table 3).

Traversing beams with force feedback
The simulation results showed that the force feedback strategy worked well, and the traversal was successful for both flimsy and stiff beams (figure 7; movie 1).The mechanical energy cost to traverse flimsy beams was 15.0 mJ, which was almost the same as that of the feedforward pushing strategy (14.8 mJ; section 3.1.1).The slight increase in energy might be attributed to its cost in pitch control to track desired pitch angle.The mechanical energy cost to traverse stiff beams was 17.4 mJ, which was much smaller than that of the feedforward pushing strategy (185.2 mJ; section 3.1.1).It also consumed less mechanical energy than the obstacle avoidance strategy (35.8 mJ; section 3.1.3),because the robot could lean against the beams between the gap and the interaction/forces torques with beams helped the robot roll.
In comparison to strategies without force sensing, our force feedback strategy had the advantage of being able to adaptively traverse beams in response to beam properties (table 4).Without force sensing, it was impossible to differentiate between beam types, leading to the use the same strategy for both stiff and flimsy beams.While the feedforward pushing strategy was effective for flimsy beams, it was challenged by stiff beams, leading to significant energy costs if the force/torque limits were not reached, or causing the robot to become stuck upon reaching that limits.On the other hand, the obstacle avoidance strategy, although working well for stiff beams, was not optimal for flimsy ones due to high energy costs.By incorporating force sensing, our force feedback strategy required minimal energy cost during the traversal of both beam types and reduced the risk of the robot getting stuck.Unlike the feedforward pushing strategy, which failed to traverse stiff beams with force/torque limits, our force feedback strategy successfully traversed them even with force/torque limits.

Traversal performance with different sensorimotor delay
The longer the sensorimotor delay was, the more the robot had pitched up when it started active control (figure 8(A)).Despite this, the robot traversed in all trials (movie 1).We found that the longer the

Discussion
Our modeling and simulation results showed that environmental force sensing can help robots determine previously unknown beam properties and inform planning and control accordingly to traverse cluttered obstacles by effectively using physical interaction.We compared our new force feedback traversal strategy with commonly used traversal strategies, feedforward pushing and obstacle avoidance, in simulation.The force feedback strategy based on force sensing helped the robot save mechanical energy cost during beam traversal, prevented it from becoming stuck in front of beams, and could prevent it from flipping over.Our modeling study demonstrated the usefulness of body contact force sensing for estimating terrain properties and allowing robots to traverse large, cluttered obstacles with enhanced performance.This novel capability to use the body to sense environmental contact forces to infer environmental physics and guide motion control will be directly useful for RHex-class robots [50,51] moving in complex 3D terrain, and it will also catalyze more diverse robot platforms to sense and use body-obstacle physical interaction (in addition to sensing leg-ground interaction [37][38][39][40][41][42]) to better move in the real world.Our lab is currently developing a robophysical model of the body-beam interaction system, instrumented with custom internal 3D force sensors and distributed external contact sensors [52].Building on the modeling insight here, we will use it to further study the principles of mechanical sensory feedback control and demonstrate advancement in a real robotic system.
In addition to robotics, our work also has implications for biology.The discoid cockroach (and other animals that use body-environment physical interaction to traverse cluttered large obstacles) uses loadsensitive mechanoreceptors distributed throughout their body [53] to sense contact.Our modeling suggested that this sensory information may help them use locomotor modes that are less costly to traverse.Such mechanical sensory feedback can complement eyes [54] and antennae [55] that sense the visual/geometric information of the environment and help them better traversing cluttered obstacles with reduced energy cost.
Our motion planning has limitations that require further examination.We neglected the effect of body inertia and assumed constant forward speed in motion planning.These approximations and assumptions simplify the system modeling and analysis, but they may compromise traversal performance.We attempted to apply optimal control while considering body inertia and allowing forward speed to vary.However, it was challenging to find the global minimum because of the non-smooth nature of contact force (for example, abrupt changes in contact force during the transition from roll to pitch mode) and intermittent collisions between the robot body and beams.
While this work simplifies obstacles by modeling them as rigid plates with torsion springs, our proposed approach should extend to a variety of other obstacles.The only prerequisite for our method is the availability of a model that establishes the relationship between contact force and physical properties of obstacles.If the robot encountered a new obstacle, we could develop a physics model tailored to it.For example, when dealing with movable rigid obstacles, considerations such as resistive force and displacement in all directions would become crucial factors in the model and could be monitored through continuous force sensing.Such an extension will be useful in search and rescue operations amidst earthquake rubble, where obstacle avoidance may not be feasible and obstacles vary in stability, rigidity, and mobility.By employing force sensing, we can formulate a traversal strategy to decide whether to lean on, avoid, or push obstacles away.
We focused on using the predictive physics model to estimate unknown obstacle physical properties.Future work should also expand its use in robot control.For example, feedback control (equation ( 7)) requires real-time reading of the robot's states and external forces.However, because there is a time delay in sensors, the actual control input is always calculated using the state and external forces from the previous moment.Therefore, when velocities or accelerations are large, the control may not be timely and may lead to unexpected collisions.In this case, the physics model can be used to predict the evolution of the state and external forces and enable more effective and timely control of robots.

Figure 1 .
Figure 1.Discoid cockroach and robophysical model traverse grass-like beam obstacles using physical interaction.(A) For flimsy beams, the discoid cockroach predominantly uses a pitch mode (blue), while for stiff beams, it predominantly transitions (yellow) to a roll mode (red)[35].(B) A minimalistic robophysical model of cockroach beam traversal[35].Adapted with permission from[35].

Figure 2 .
Figure 2. Multi-body dynamic simulation of minimalistic robot traversing beam obstacles.(A) The simulation robot body has five degrees of freedom, fore-aft, vertical, lateral, roll (red arrow), and pitch (blue arrow), which are controlled by motors.The roll axis (red dashed line) is an axis through the rotation center of the robot and parallel to the X-axis of the world frame.The pitch axis (blue dashed line) is the y-axis in the body frame.An additional weight makes it bottom-heavy.Each beam is connected to the ground via a torsion spring and can deflect in the vertical X−Z plane.The world frame is X-Y-Z.Body frame is x-y-z, which is attached to the robot.(B), (C) Snapshots of simulation and physical robot experiments from our previous study [35].(B) Without vertical oscillations, both simulation and physical robots are stuck in the pitch mode.(C) With sufficient vertical oscillations, both roll into the beam gap and traverse using the roll mode.

Figure 3 .
Figure 3. Physics model of robot-beam interaction.(A) Cockroach-inspired robot is modeled as an ellipsoidal rigid body.FX is propulsive force.τ 1 and τ 2 are roll and pitch torques along the roll and pitch axes in figure 2. The center of mass (C) is lower than the geometric center (O, the intersecting point of roll and pitch axes) by hc.(B) Each beam is modeled as a rigid plate with a torsion spring with stiffness k. (C) ⃗ F1 and ⃗ F2 are theoretical contact forces with each beam calculated from the model.(D) Left: in the pitch mode, the potential energy is at a maximum when the bottom of the robot reaches the top of the beam.Right: in the roll mode, the potential energy is at a maximum when the robot rolls until there is no beam contact.

Figure 4 .
Figure 4. Motion planning.(A) 3D potential energy landscape of the robot-beam interaction system.(B) Schematic of the weighted directed graph based on potential energy landscape.Circles are states.Values in circles show the potential energy (in Joules) at each state, while values on directed edges show the predicted mechanical energy cost of the movement along the edge.(C) The path from the initial to the target state in the weighted directed graph, with the minimal predicted mechanical energy cost.

Figure 5 .
Figure 5. Representative snapshots of simulation robot traversing stiff and flimsy beams using different strategies.(A) Feedforward pushing strategy with unlimited propulsive force.The robot traversed both stiff and flimsy beams.(B) Feedforward pushing strategy with limited propulsive force.The robot was stuck in front of stiff beams but pushed over flimsy beams.(C) Obstacle avoidance strategy with force/torque limits.The robot successfully traversed both stiff and flimsy beams without contact.

Figure 6 . 2 )
Figure 6.Mean of relative estimation error ( ek1+ek2 2 ) as a function of oscillation amplitude.Five trials were performed at each oscillation amplitude.The blue curve represents the mean, and the shaded area indicates the standard deviation (±s.d.).

Figure 7 .
Figure 7. Snapshots of the robot traversing stiff and flimsy beams using force feedback.

Figure 8 .
Figure 8. Influence of sensorimotor delay.(A) Snapshots showing system configuration when the control started after different sensorimotor delays.(B) Mechanical energy cost from forwarding and rotation (both roll and pitch) actuators (FX and TR cost) as a function of sensorimotor delay (one trial each).

Table 1 .
Error in beam stiffness estimation (in N m rad −1 ) for varying oscillation frequencies (five trials each).

Table 2 .
Error in beam stiffness estimation (in N m rad −1 ) for varying oscillation amplitudes (five trials each).

Table 3 .
Error in beam stiffness estimation for different sensing times (five trials each, f = 2 Hz, A l = 1 mm).

Table 4 .
Mechanical energy cost of different strategies for traversing stiff and flimsy beams (one trial each).
sensorimotor delay, the more energy cost from forwarding and rotation (both roll and pitch) actuators (figure8(B)).Thus, the faster the robot could estimate the properties of obstacles and respond after encountering obstacles, the more energy it would save.