Frictional Locomotion of a Radially Symmetric Tripedal Robot

This study seeks to provide physical insight into the friction-driven crawling locomotion of systems with radially symmetric bodies. Laboratory experiments with a tripedal robot show that both translation and rotation can be achieved with just three independently actuated rigid limbs, i.e., with 3 degrees-of-freedom. These observations are rationalized using a simple mathematical model, which assumes that the friction at each limb is linearly proportional to the normal force at the contact point, and opposes the direction of motion. This dynamic model reproduces experimental observations across an extensive parametric sweep involving sinusoidal rotation of the limbs with varying amplitudes and phase shifts. Model predictions highlight the role played by time-varying normal forces at the contact points. These predictions are confirmed using embedded force transducers in the limbs. We present a further simplified analysis explaining that a geometric nonlinearity is induced in the dynamics from the radial symmetry and that this nonlinearity is essential to the generation of pure translation. We also show that this nonlinearity can be amplified by a cyclic time-varying limb length variation. These results provide a framework for further study of radially symmetric movers.


Introduction
Previous locomotion studies have primarily focused on bilaterally symmetric body arrangements. This is partly because the majority of land-based animals are bilaterally symmetric. However, radially symmetric body arrangements are common in aquatic environments, e.g. Cnidaria (swimming medusae, or jellyfish; Cameron and Fankboner (1989)) and Echinodermata (crawling organisms such as sea stars, brittle stars, and sea urchins; Manuel (2009); Astley (2012)). We are interested in exploring how these predominantly aquatic body arrangements can be used for amphibious or terrestrial locomotion. As a starting point in this endeavor, this study analyzes friction-driven locomotion of tripedal geometries involving three contact points-the minimum number required for stable motion with radially symmetric bodies-via laboratory experiments and the development of simplified mathematical models. The focus on a tripedal geometry allows us to minimize the complexity of the robot used in the experiments and provide simplified analyses for the nonlinear dynamics we observe.
Radially symmetric morphologies offer some benefits over bilaterally symmetric arrangements. First, radially symmetric movers can use any limb to define a lead direction Cole (1913); Astley (2012). For example, pentaradial brittle stars are shown to exhibit bilaterally symmetric motion strategies in which one of the five limbs remains inactive Astley (2012). However, because they are radially symmetric, the inactive limb can be reassigned to change direction immediately Astley (2012). In addition to having no front bias, in the case of radially symmetric body plans with more than 3 limbs, the mover is robust to limb damage due to redundancy Kano et al. (2017Kano et al. ( , 2019; Sköld and Rosenberg (1996); Carnevali (2006); Matsuzaka et al. (2017). This advantage may be of use in remote environments where loss of motor function or general damage to a limb does not render the mover immobile. Finally, for locomotion in the presence of strong external forces (e.g., aerodynamic or hydrodynamic disturbances), having multiple contact points enables more robust adhesion. This can reduce risk of surface dislodgement and enhance the net adhesive force through distributive attachment Arshavskii et al. (1976).
Several recent efforts have considered the development of radially symmetric robots. Our previous work describing the development of the tripedal test system used in this study showed that the radially symmetric arrangement is controllable using an omnidirectional gait map Hermes et al. (2021). This prior effort also highlighted an advantage of being able to define arbitrary lead limbs for radially symmetric robots: sharp path curves can easily be followed. The Robotics and Mechanisms Laboratory at UCLA has produced many radially symmetric systems: the hexapod robots SiLVIA, HEX, and MARS Showalter et al. (2008); a radially symmetric quadruped, ALPHRED Hooks et al. (2020); and two tripod robots, THALer Webb et al. (2015) and STriDER Heaston et al. (2007). In general, these robots have multi-actuated legs for walking and interacting with their environments. The tripedal STriDER and THALer systems have also demonstrated application of a novel swinging gait allowing for translation in three directions. The Ishikawa lab group Oki et al. (2015); Ishikawa et al. (2018); Masuda and Ishikawa (2017) have also demonstrated controlled crawling using a robot called Martian with similar geometry to STriDER. Unlike STriDER, this robot uses small vertical stepping oscillations coupled with joint rotation driven by a decentralized control scheme using normal force feedback. Locomotion studies for this robot showed natural synchronization to occur from a large set of initial conditions. Bevly et al. (2000) have demonstrated vertical climbing using a radially symmetric tri-arm robot. Several other groups have developed hex-symmetric terrestrial Faigl andČížek (2019); Žák et al. (2021); Liu et al. (2018); Bjelonic et al. (2018) or aerial Ryll et al. (2017) robots.
Other research efforts have considered the development of robots with close similarity to pentaradial brittle stars. For example, Watanabe et al. (2012) developed a decentralized control scheme to orchestrate limb motion for a pentaradial system. Kano et al. (2012) demonstrated omnidirectional locomotive capabilities using this control scheme on a physical ophiuroid robot. More recently, Kano et al. (2017) demonstrated robust locomotion capabilities in a pentaradial robot by systematically removing appendages and showing the adapted locomotion strategy. This effort built on biological observations that involved removing arms from live brittle stars and observing the resulting locomotion patterns. Lal et al. (2007Lal et al. ( , 2008 developed motion control algorithms for a pentaradial robot with multisegmented arms using genetic algorithms. This work showed that a nonintuitive complex writhing motion generated the largest net translation. Given the triradial geometry considered in this study, a related field of research involves controlling the motion of trident robots, which are tripedal robots with wheels at the end effector instead of simple frictional contacts Ishikawa et al. (2009a, b); Jakubiak et al. (2010);Yamaguchi (2012). In general, these studies use nonholonomic kinematic relations to obtain optimal control strategies. Ishikawa et al. Ishikawa et al. (2009a) have also demonstrated the ability to track sharp path trajectories with a geometric mechanics algorithm.
The efforts described in the preceding paragraphs highlight the broad applicability of radially symmetric robots, and emphasize the comparative advantages to systems with bilateral symmetry. The present work seeks to build on these prior efforts and provide insight into the nonlinear dynamics underlying friction-driven sliding (or crawling) locomotion of systems exhibiting triradial symmetry. For this purpose, we have developed an idealized robot with 3 rigid single DOF limbs. We have characterized the motion of this robot in laboratory experiments comprising a systematic parametric sweep involving open-loop sinusoidal actuation of the limbs. We have also developed a simplified physics-based mathematical model that reproduces the experimental observations, and provides insight into the unique locomotive capabilities of this arrangement. Through this effort, we have identified and characterized actuation patterns that lead to pure rotation and pure translation.
The remainder of this paper is structured as follows. Robot design and laboratory experiments are described in Sect. 2. The physics-based model is developed in Sect. 2.3. Experimental results and model predictions are compared in Sect. 3. Actuation patterns that lead to pure rotational and translational motion are also evaluated in this section. We then seek an explanation for answering why the symmetric gait sequence generates net translation in Sect. 4. Concluding remarks are presented in Sect. 5.

Experimental Approach and Physics-Based Model
In this section, we provide an overview of the robot design (2.1), our mathematical model (2.3), and normal force validation (2.4).

Robot Design and Characterization
A top and isoview of the simple tripedal robot used in this study is shown in section 1. The tripedal robot comprises 7.5cm long 3D-printed polylactic acid (PLA) arms connected to an stereolithography (SLA) resin motor support structure with an effective radius of 5cm. We used Hitec HS-5646 digital servo motors to drive the arms, which were powered by an external power supply and controlled by an Arduino Uno. The total mass of the structure was measured to be 888g. We also embedded limbs with 5 kg bending-beam load cells (Chenbo CZL635) with an HX-711 amplifier module to measure normal forces with a noise profile of roughly ±1.5%FS (approximately ±0.7N).
We varied the friction at the contact points using polymer laminate and 600-grit sandpaper surfaces. We used paper as the test platform surface. The friction coefficients measured were μ = 0.33±0.012 for the polymer-paper contact and μ = 0.85±0.043 for the sandpaper-paper contact. As described in Hermes et al. (2021), these values were obtained by releasing a known mass (with appropriate contacts) from rest down a slope and measuring the velocity achieved by the mass after a set distance. The difference between the initial gravitational potential energy and the final kinetic energy of the sliding mass was attributed to friction, and friction coefficients were estimated from this frictional loss. The sliding masses were randomly oriented and produced small deviations in the measured friction coefficients. From this, we concluded that the friction coefficient is nearly isotropic. As discussed below, we ignored high-order friction effects for modeling simplicity. To ensure that the observed robot motion was purely due to limb actuation, we ensured that the dynamic impact of surface imperfections and the power/control tether were minimized. All experiments were carried out on a smooth, level optimal support table with tilt less than 0.2 • . The tether cable was supported externally to minimize tension.

Motion Experiments and Video Tracking
To characterize the dynamic behavior of the tripedal robot, we pursued a systematic experimental evaluation in which we measured the robot motion resulting from prescribed sinusoidal actuation of each limb. Specifically, the angular position of the i−th limb (see Fig. 2) was prescribed as where i = 1, 2, 3. Here, ξ represents rotation amplitude, f is the oscillation frequency, and ψ is a temporal phase shift in actuation between the limbs. We refer to ξ as sweep and ψ as phase for the remainder of this paper. The following ranges were tested for each of the three actuation parameters: The ranges for f and ξ were constrained by servomotor limitations and stability considerations. The range of test values for ψ was chosen to span the configuration space. Table 1 shows the specific parameter combinations tested. In addition to this sinusoidal parametric sweep, we also tested actuation patterns resembling locomotion strategies observed for brittle stars named rowing and reverse rowing Astley (2012). Both of these translation modes involve one unpaired or inactive limb. Rowing and reverse rowing are distinguished by the orientation of this inactive limb and are graphically described in Fig. 9. In the case of rowing, the inactive limb leads the body forward, whereas for reverse rowing, the inactive limb passively trails the body. For both locomotor   Astley (2012). Thus, despite the radial body plan, brittle stars effectively employ coordinated, bilaterally symmetric locomotion; radial symmetry is broken by the presence of an inactive limb. The robot was tracked at 30 frames per second from above using a camera mounted on an elevated platform. Because the largest test frequency was 1.2 Hz, the selected frame rate resolves motion completely. The center of mass for the robot in the horizontal plane, x = [x, y], was determined by applying a MATLAB image-filter script to the recorded images. A reference measurement object was used for pixel-to-centimeter calibration in the motion plane. The relative rotation of the robot, θ , was tracked by finding the maximum correlation between the preceding image, rotated in the range [−3 • , 3 • ], and the frame being considered. The starting position of the robot for each experiment was defined as x = [0, 0] and θ = 0 • .

Mathematical Model
We developed a simplified dynamic model for the tripedal robot based on momentum conservation laws with nonlinear friction forces. In our model, we assume rotational and translational intertias for the limbs to be negligible compared to those for the central body. Friction forces are modeled as point forces at the tips of the limbs.
The friction force is simplified in our model to be the signum function of the velocity. In other words, there are no stick-slip effects. Olsson et al. (1998); Pennestrì et al. (2016). We recognize that this friction model constitutes a significant simplification. However, below we show that this simplified model is able to reasonably reproduce experimental observations. The friction force, F i = [F i,x , F i,y ], acting at the contact point at the end of each limb in the x − y plane, r i = [r i,x , r i,y ], is modeled as Here, N i is the normal force at contact point i, μ is the measured kinetic friction coefficient, is the velocity of the contact point at the end of each limb. Note that r i andṙ i can be expressed in terms of the robot state vector, z = [x,ẋ, θ,θ ], the actuation angles and rotation rates, φ i andφ i , and the geometric constants, R and l, using simple trigonometric relations (see Fig. 2). To avoid numerical issues for cases in which the velocity at the contact point is zero, we assume F i = 0 for |P r i | 1, or that there is no translational friction force for infinitesimally small velocities. Normal forces are computed by combining a vertical force balance (N 1 + N 2 + N 3 = Mg) with torque balances about the robot center of mass: The second and third lines in the equation above ensure that there are no net torques about the robot center-of-mass due to the normal forces. Under the assumptions stated above, conservation laws for linear momentum in the horizontal plane can be expressed as and respectively. The conservation law for angular momentum can be expressed compactly as where J represents the second mass moment of inertia, which is estimated from the mass distribution of the central support structure.
The system of ordinary differential equations shown in (4-6) is solved numerically to yield predictions for robot translation (x) and rotation (θ ) using the MATLAB ode45 algorithm for an adaptive time-step 4th-order Runge-Kutta solver. Recall that the location of the contact points relative to the robot center, r i −x, depends on the limb angles, φ i (t) (Fig. 2). Therefore, the prescribed actuation angles φ i (t) appear directly The oscillation frequency and sweep angle are set to f = 0.5 Hz and ξ = 30 • for the high-friction sandpaper-paper contact. Shaded red regions represent measurement variability. This is defined as one standard deviation for the phase averaging procedure in all three conservation laws via the friction terms dependent on r i andṙ i . Also, keep in mind that these model simulations involve no tuning parameters. All geometric and dynamic variables appearing in the governing equations (e.g., M, J , l, R, μ) are obtained from independent measurements.

Normal Force Simulation Predictions and Measurements
As noted earlier, we embedded three load cells into the robot limbs to measure normal forces. In this section we compare the measured normal forces against predictions made using eq. (3). Measured normal forces for f = 0.5Hz, ξ = 30 • , and varying phase differences ψ = 60 • , ψ = 90 • , and ψ = 120 • are plotted with simulation data in Fig. 3. Note that the measured data are phase averaged over 5 oscillation cycles.
Though the measurements are qualitatively similar to the predictions, there are some consistent quantitative discrepancies. First, the amplitude of force oscillation for all limbs is larger for the measurements compared to the simulations. In addition, for ψ = 120 • , the normal force profiles of the legs are not equal with a phase shift. Unlike the simulations, the measured normal forces differ in magnitude. These measurements suggest that the robot weight may not be evenly distributed, which is a source of potential modeling error. This error could, in part, explain some of the discrepancies in trajectories between simulation data and experiment results observed in the next section. Nevertheless, the qualitative agreement between model predictions and measured normal forces is promising, especially keeping in mind the expected ±0.7N noise floor in the load cell measurements, which is comparable to the discrepancy between measurements and predictions.  Table 1. Comparison between experimental tracking results and model predictions for these cases are shown in Sects. 3.1-3.3. Simulation predictions were interpolated at a rate of 30 Hz, consistent with the experimental imaging results. The total time interval for all the results shown in Figs. 4-6 below was 30 s. Sections 3.4-3.5 build upon these findings to consider actuation patterns that result in pure rotation and translation.

Frequency Variation
Fig. 4 compares model predictions for robot position with tracking results from experiments with varying oscillation frequency, f . The sweep angle and relative phase shift were maintained constant at ξ = 30 • and ψ = 30 • for these experiments. Model predictions and experimental results both show the presence of small-scale oscillations in robot center-of-mass as well as a larger-scale orbit or revolution. For the low-friction cases shown in (a)-(d), an increase in frequency led to a contraction in the local oscillation length in the experiments, i.e., the small-scale oscillations clustered closer together. For the high-friction cases shown in (e)-(h), increasing frequency did not significantly change the local oscillation scale or the arc radius of the large-scale revolution. Instead, increasing frequency simply led to an increase in the total distance traveled by the robot, i.e., the arc length traversed by the robot increased.
For the high-friction cases, model predictions (black lines) are in good qualitative and quantitative agreement with experimental observations (red lines). For the lowfriction case, the model predictions are in qualitative agreement with the experimental Panels (a-d) correspond to the low-friction tape-paper contact while panels (e-h) correspond to the high-friction sandpaper-paper contact. Tracking data from the experiments are shown as thin red lines and the simulation data are shown as a thick black lines results, though the clustering effect is not reproduced as clearly; only a minor reduction in the oscillation scale is observed from panel (a) to panel (d). Note that the discrepancy between model predictions and experiments is greatest for the lowest-frequency cases considered in (a) and (e). This discrepancy could be attributed to the highly-simplified contact model considered here, which does not distinguish between static and kinetic friction. Transitions between static and kinetic friction (i.e., stick-slip dynamics) are expected to be more important at lower actuation frequencies.

Sweep Variation
Next, we consider robot motion for cases in which the sweep values were varied systematically while the oscillation frequency and phase shift were kept constant at f = 1.0 Hz and ψ = 30 • . Experimental observations and model predictions of robot trajectory for these cases are shown in Fig. 5. For both the low-and high-friction tests, an increase in sweep angle led to an increase in the local oscillation length. For the high-friction cases, the arc radius for the large-scale revolution also decreased slightly with increasing sweep angles. Together, this increase in oscillation length and reduction in arc radius resulted in a lower revolution period. In other words, for the high-friction cases, the robot completed a large scale revolution faster as the sweep angle increased. This is clear from a comparison between the high-sweep case shown in Fig. 5(h), which shows a complete revolution, and the lower-sweep cases shown in Figs. 5e-g, which do not show a complete revolution. The low-friction cases do not show a clear reduction in arc radius, though the total distance traveled by the robot does increase monotonically with increasing sweep angles. Once again, the mathematical model generates predictions in good quantitative agreement with the experimental observations for the high-friction cases. For the lowfriction contact, the model qualitatively reproduces the observed increase in the local oscillation scale and distance traveled by the robot with increasing sweep. However, there are deviations between the observed trajectories and model predictions.

Phase Variation
Finally, Fig. 6 shows the effect of varying phase difference, ψ, on robot motion at constant oscillation frequency, f = 1.0 Hz, and sweep angle, ξ = 30 • . For both the low-friction and high-friction cases, an increase in the phase shift led to a decrease in the arc radius for the large-scale revolution. This was accompanied by a marked change in the character or shape of the small-scale oscillation patterns. Specifically, the oscillations became sharper and less rounded with increasing phase separation between limbs. This is in contrast to the results shown in Fig. 5, in which the arc radius decreased slightly with increasing sweep, but the rounded profile of the oscillations remained. As with the sweep and frequency variation experiments, there is very good quantitative agreement between model predictions and experimental trajectories for the high-friction cases, but less so for the low-friction cases.

Pure Rotation
The observed decrease in arc radius with increasing phase shift observed in Fig. 6 inspired an additional experiment in which the phase shift was set to ψ = 120 • . Given the triradial symmetry of the robot, for ψ = 120 • , the sinusoidal actuation described by (1) leads to a constant phase difference between limbs i = (1, 2) and i = (2, 3) but not i = (3, 1). For ψ = 120 • , the phase difference is constant across all limbs. Figure 7 shows the time-variation in the x-location of the robot center of mass (a) and rotation angle θ (b) for sinusoidal actuation with ψ = 120 • . The observed trajectory shows minimal translation (x ≈ 0 m) and a monotonic increase in the rotation angle, i.e., the robot effectively exhibits pure rotation for ψ = 120 • . Model predictions for translation in the x-direction are similar to the experimental observations. However, they differ in predicting both the average rate and nature of the rotation. The model predicts a pronounced acceleration-deceleration cycle in rotation rate at the actuation frequency ( f = 1.0 Hz), whereas the experiment results show a relatively constant rotation rate.
To provide further insight into the observed rotational -or pivoting -motion, additional simulations were pursued for ψ = [20, 30, ..., 350] • , with the sweep angle and oscillated frequency fixed at ξ = 30 • and f = 1.0 Hz. The resulting largescale revolution was characterized by fitting a circle with radius r and center (a, b) in the horizontal plane to the observed trajectories. Specifically, a quasi-Newton gradient optimization method (with tolerance 10 −6 ) was used in MATLAB to minimize the error function for all data points (x j , y j ) in the predicted trajectories. The normalized curvature obtained from this procedure, l/r , is plotted in Fig. 8 as a function of the phase different ψ. The predictive curvature is maximum (i.e., fitted radius r is minimum) at ψ = 120 • and 240 • . As noted earlier, given the triradial symmetry of the robot, ψ = 120 • ensures a constant phase difference across all three limbs. This is also true for ψ = 240 • . Interestingly, the model also predicts zero curvature (i.e., r → ∞) for ψ = 180 • . This particular actuation essentially corresponds to a rowing motion with an oscillating unpaired limb. As we discuss in the following section, this actuation results in translational motion with minimal curvature. Model predictions shown in Fig. 8 indicate that the large-scale revolution is driven by the inconsistent phase differences resulting from ψ = 120 • (or 240 • ). The physical mechanism giving rise to this effect is illustrated by Fig. 3, which shows the time-varying normal forces acting at each contact point for phase shifts ψ = [60, 90, 120] • . As expected, in all cases, the normal forces show periodic variation at the oscillation frequency. However, the cases with ψ = 60 • and 90 • are characterized by the presence of a single differentiated limb that exhibits a higher-amplitude oscillation (see black lines in Figs. 3a, b). This differentiated limb corresponds to i = 2 in (1), which maintains the prescribed phase shift, ψ relative to the other two limbs. Recall that the phase shift between the remaining two limbs, i = 1 and i = 3, is not ψ except for cases in which ψ = 120 • (or a multiple thereof). These differences in relative phase give rise to differences in the frictional forces acting at each limb, which drive the large-scale orbital motion exhibited by the robot. For ψ = 120 • , there is no differentiated limb. As shown in Fig. 3c, the normal force waveforms are equal in amplitude and exhibit a constant phase difference relative to one another. In this case, the robot effectively rotates in place.

Pure Translation: Brittle Star Inspired Rowing & Reverse Rowing
In addition to characterizing the motion resulting from the sinusoidal actuation in (1), we were also inspired to reproduce brittle star locomotion patterns observed by Astley Astley (2012), which we interpret as the rowing and reverse-rowing motions shown in Fig. 9. Exact replication of the complex kinematics exhibited by brittle stars is impossible with the idealized rigid robot. However, the rowing and reverse-rowing motions can be simulated qualitatively by making one of the robot limbs inactive, and by altering the neutral position and motion of the remaining two active limbs. To simulate rowing, the neutral position of the limbs was unchanged from the natural state in which the limbs are distributed evenly around the unit circle, i.e., separated by 120 • from each other. To simulate reverse-rowing, the neutral position of the active limbs was shifted by 30 • towards the inactive limb. As shown in Fig. 9, this meant that the active limbs were each separated from the inactive limb by 90 • in the neutral position, and from one another by 180 • . For both rowing and reverse rowing, the active limbs were actuated in anti-phase, i.e., with a phase difference of 180 • . Figure 9 provides a comparison of experimental tracking results and model predictions for the simulated rowing and reverse-rowing motions with the following actuation parameters: oscillation frequency f = 1.0 Hz and sweep angle ξ = 30 • . The rowing results shown in panels (b) and (d) show sustained translation in the positive x-direction. In other words, the robot moves in the direction of the inactive limb, which is consistent with biological observations. Model predictions are also in good agreement with the experimental tracking data. In contrast to the rowing locomotion results presented in Fig. 9b and d, the attempt at simulating reverse-rowing was unsuccessful. As shown in panels (a) and (c), the robot effectively oscillated back and forth in place, with no significant net motion relative to the starting point at x = 0. Model predictions also show similar behavior. Thus, simply altering the neutral position of the active limbs does not reproduce the reverse-rowing locomotion observed in nature. This issue is discussed in greater in the following section.
To evaluate the effect of varying actuation parameters on net translation speed for the rowing configuration, we pursued additional experiments and model simulations with sweep angles ξ = [20,30,40,50] • and oscillation frequency f = 1.0 Hz. In other words, the oscillation amplitude for the active limbs was varied. Figure 10 shows the net translation of the robot center-of-mass as a function of time, x(t). In general,   To provide further insight into the observed reduction in net translation speed for the high-friction sandpaper-paper contact, model predictions were carried out with friction coefficients μ ∈ [0.1, 0.9] and sweep angles ξ = [20,30,40,50] • . The average translation speeds predicted by the model are shown in Fig. 11. Interestingly, for each value of ξ , there is an optimal friction coefficient that maximizes translation speed. Beyond this optimum, translation speed decreases monotonically with increasing friction coefficient This is consistent with the experimental observations in Fig. 10, which showed a decrease in average translation speed for the high-friction sandpaper-paper contact. Model predictions also suggest that an increase in the sweep angle increases the maximum translation speed and shifts the optimum friction coefficient higher. Physically, the presence of an optimum could be explained by considering the following limiting scenarios. For μ = 0, the active limbs would simply slide back and forth on the substrate without generating any net friction forces (or translation). However, for very high friction coefficients, the inactive limb may effectively act as an anchor that slows the robot down.

Analysis of Translation
In this section, we pursue a simplified analysis to provide further insight into rowing and reverse rowing locomotion. Specifically, we seek to provide an explanation for why the attempt at replicating reverse rowing was unsuccessful with the rigid-limb tripedal robot used in this study. Without loss of generality, the representation shown in Fig. 12 can be used to describe the kinematics of one active limb for the rowing and reverse rowing modes considered in Sect. 3.5. The horizontal and vertical components of velocity for the active limb,ṙ = [ṙ x ,ṙ y ], can be written aṡ whereẋ = [ẋ,ẏ] represents the velocity for the robot center of mass, α ∈ (0, π) represents the constant offset angle for the active limb and φ = ξ sin(ωt) represents the angular actuation with sweep angle ξ and frequency ω = 2π f . Per the model shown in eq. (2), the friction force generated by the active limb in the direction of motion is expected to be proportional to Given the symmetric actuation of the active limbs in rowing and reverse rowing locomotion, the other active limb would generate an identical frictional force in the direction of motion. The above analysis assumes that the normal force is constant. To test the validity of this assumption, we performed simulations with a constant normal force and found the character of the resulting motion patterns to not change. Modeling the time-varying normal forces is important for the accuracy of the numerical simulations. However, the same trends can be captured with a constant normal force approximation.
For there to be net motion in the horizontal direction, the cycle-averaged value of the friction force must be non-zero,F x = 0. To evaluate the conditions in which this occurs analytically, we redefine the offeset angle as α = π/2 + and consider the limit in which the actuation angle is small, |φ| 1. Under this assumption, we have such that To make further analytical progress, we assume that the horizontal motion of the center of mass can be expressed as: in which A and B represent a scaling factor and phase shift in the motion of the center of mass relative to the leading harmonic in the actuation velocity, φ sin α = ωξ sin α cos(ωt). Note that the above expression neglects the (small) cycle-averaged motion of the center of mass as well as any motion at higher harmonics. Substituting the assumed expression forẋ into the simplified friction model and factoring out ω sin α yields (14) The above expression indicates that the following three conditions are necessary for the generation of a non-zero cycle-averaged friction force: (1) = 0, (2) A = 0, and (3) B = nπ with n ∈ Z. This can be shown by considering each of the three conditions separately.
(1) If = 0, the friction force expression simplifies to: or for which F x = 0 since the signum function for waveforms comprising single Fourier harmonics is positive and negative for the same duration of an oscillation cycle.
It can be shown formally (or via simple numerical tests) that the expression above again yields zero cycle-averaged force, F x = 0.
Condition (1) explains why reverse-rowing, which is characterized by an offset angle α = π/2 or = 0 in the present work, is unsuccessful. Physically, the geometric nonlinearity in limb motion that arises for α = π/2 and gives rise to a second harmonic in horizontal limb velocity,ṙ x ∝ sin(2ωt), is essential to the generation of a cycleaveraged friction force. The results presented in Sect. 4.2 show how this geometric nonlinearity can be reintroduced via a change in limb length over an actuation cycle.
Conditions (2) and (3) indicate that there is no cycle-averaged force if the center of mass is held in place (A = 0) or if the motion of the robot center of mass is exactly in phase or antiphase with the motion of the active limbs. Physically, for dynamics governed by Eq. (4), these conditions are expected to be satisfied naturally. In other words, we expect non-zero motion of the robot center of mass that is out of phase with the motion of the limbs. Indeed, as shown in Fig. 13, the leading harmonic in robot x-direction body velocity is described for simulations by scaling factor A ≈ 0.30 and phase shift with x-direction limb velocity, B ≈ 0.08 [rad]. For experiments, the scaling factor is A ≈ 0.57 and the phase shift is B ≈ 0.23 [rad]. Note that the experimental velocities were estimated using a first-order derivative of position data obtained from frame-by-frame tracking.

A Successful Reverse Rowing Gait
The results obtained in the previous section showed that reverse-rowing is not possible with the robot developed and tested in this study due to condition (2). This suggests that, in general, rowing locomotion is likely a more effective strategy for radially symmetric crawling robots with fixed limbs. To test if changing limb length could indeed enable crawling with the reverse-rowing mode, we pursued additional simulations with the mathematical model from Sect. 2.3. Note that the results described in the previous sections indicate that this model generates very useful qualitative and quantitative predictions, even if perfect agreement with the experiments is not observed in all cases.
To mimic the reverse-rowing sequence exhibited by brittle stars where arm length varies over a cycle Kano et al. (2017), the length of the active limbs in our model was prescribed to vary as follows over an oscillation cycle: The parameter l s represents the minimum limb length during the cycle, and l o = 0.075 m is the maximum limb length. As shown in Fig. 14, this reproduces a breast stroke- The results shown in Fig. 15 confirm that reverse rowing is possible with varying limb length. Further, the translation speed increases monotonically with decreasing l s , i.e., as the relative change in limb length increases. The model also indicates that translation speed for the modified reverse-rowing mode increases initially with increasing friction coefficient, before approaching an asymptotic limit that depends on l s . This is unlike the results obtained for rowing locomotion, which predict that there is an optimal friction coefficient that maximizes translation speed. Further increases in μ beyond this optimum lead to a decrease in translation speed. These observations suggest that the modified reverse-rowing locomotor mode may be less sensitive to substrate type. In general, model predictions presented in this section indicate that the ability to vary limb length could lead to more robust locomotion capabilities. For robotic systems, variable limb length could be achieved either via the use of soft appendages or by incorporating additional degrees of freedom into rigid limbs, such as in Oki et al. (2015).

Non-sinusoidal Gait Optimization
This study largely focuses on the analysis of sinusoidal gait patterns. However, there may exist non-sinusoidal gaits that result in faster translation. As a preliminary evaluation of this effect, we revisited the rowing strategy described in §3.5 and pursued parametric optimization of the cyclic actuation of the limbs, i.e., the angular position φ(t) in (1).
For this optimization, we used a three degree-of-freedom continuous spline that is constrained to be zero at the beginning and end of the cycle, and also constrained to cross the zero axis only once. As a starting point, a Piece-wise Cubic Hermite Interpolating Polynomial (PCHIP) spline was used. We selected PCHIP because it allows for sharper curvatures compared to sinusoidal functions and a more diverse set of allowable functions. Limb oscillation was parameterized using three variables, p = {p 1 , p 2 , p 3 }, defined as: where T = 2π/ f . The variables p 1 , p 2 , and p 3 represent the maximum angular position, the minimum angular position, and the time corresponding to the zero crossing in the oscillation cycle, respectively. The objective function, x f (p), for the optimization was the final x−position of the robot after 5 simulation cycles (i.e., 5 s). We selected 5 cycles to eliminate transient effects, while minimizing simulation time. To maintain consistency with the earlier rowing results, the sweep angle was set to ξ = 30 • and the oscillation frequency was fixed at f = 1.0 Hz.
Isocontours of x f are shown in Fig. 16a. The highest performing parameter combination, p 1 = 0.981, p 2 = −0.974 and p3 = 0.478, resulted in nearly sinusoidal gaits. Based on this observation, we also tested similarly parameterized piecewise-sinusoid curves instead of PCHIP curves, and found the resulting x f -map to be similar with a significant computational time reduction (see Fig. 16b). Thus, for additional gait design and optimization, we used the piecewise-sinusoid curve. We applied an adaptive-step gradient search algorithm. The resulting functional was multimodal near the optimum and the optimal parameter sets were found to be near p ≈ {1, −1, 0.456}. In other words, the optimal PCHIP and piecewise-sinusoid gaits are both close to sinusoidal (see Fig. 17a). However, by moving the crossover point earlier in the cycle, net translation increases slightly (∼ 5%) for the piecewise-sinusoidal gait compared to a purely sinusoidal gait after 5 cycles (see Fig. 17b). This preliminary evaluation of optimal gaits suggests that the study of harmonic sinusoidal gaits remains a useful starting point for future studies.
Interestingly, the optimization results in Fig. 16 also show that it is possible to get backwards translation (x f < 0) by introducing temporal asymmetry in the gaits ( p 3 > 0.5). This is consistent with the findings presented in §4.2, which indicate that kinematic asymmetry is required to generate a successful reverse rowing gait.

Conclusion
This study demonstrates that several locomotion strategies are possible for robotic systems with radially-symmetric bodies. Sequential actuation of all limbs can give rise to locomotion along circular orbits with controllable radius. Periodic actuation with a constant phase shift across all limbs ( ψ = 120 • for a tripedal robot) leads to rotation in-place. Pure translation is achieved through the presence of an inactive limb that breaks radial symmetry. The ability to move along prescribed curves, rotate in place, and travel in a straight line establishes a framework for future motion planning efforts for radially-symmetric crawling robots. Using the models developed in this paper, we have successfully demonstrated path following for this robot using image-based feedback control Hermes et al. (2021). The primary advantage of radial symmetries highlighted in this work is the omni-directional translation ability. Consequently, robots that require instantaneous omni-directional actuation may benefit most from applying the ideas presented here.
For systems with rigid limbs and sinusoidal actuation, such as the robot developed here, we found only the so-called rowing locomotor mode observed in nature to be effective at locomotion (we cannot comment on feasibility with non-sinusoidal gaits). For this mode, the hind limbs are actuated in anti-phase with one another and the robot moves in the direction of the inactive forelimb. In Sect. 4, we showed via further simplifications to the dynamic model that net translation is only present if certain conditions are met. These conditions establish that a phase shift between the robot center velocity and the limb velocity is generated due to inertial effects, and that this phase shift creates a positive net force over a cycle. Reverse-rowing was only feasible for more complicated systems capable of varying limb length over an actuation cycle. Note that the ability to vary limb length ultimately allows biological and engineered systems to vary the position of the contact points relative to the center of mass. This allows for finer control over the normal and, hence, friction forces acting at each contact point.
The ability to modulate frictional forces can play an important role in crawling locomotion. Indeed, previous studies show that snakes try to maximize normal forces over low-friction surfaces through several different mechanisms; this includes lifting parts of their body off the substrate to enhance normal forces Hu et al. (2009);Marvi and Hu (2012). Predictions made using the simple Coulomb friction model indicate the presence of an optimal sweep-dependent friction coefficient that maximizes translation speed for brittle star inspired rowing locomotion. The optimal friction coefficient values were μ < 0.3 for the range of sweep angles tested here. For higher friction coefficients, translation speed decreased. However, model predictions also suggest that this deterioration of performance can be mitigated by switching to a reverserowing strategy, in which limb length varies over an actuation cycle. These differences in crawling efficacy may explain the prevalence of both rowing and reverse-rowing locomotion in nature Astley (2012). The relative performance of these rowing and reverse-rowing locomotor modes in terms of speed and energetic efficiency remains to be evaluated.
Finally, it must be emphasized that the simple model developed here is very much a starting point for further research. More sophisticated frictional models that account for stick-slip dynamics Heslot et al. (1994), or analytical models that take advantage of scale separation between the timescale associated with the fast periodic actuation and the slow macro-scale motion (i.e., translation, large-scale revolution) might provide significant further insight into crawling dynamics. Models grounded in geometric mechanics Kelly and Murray (1995); Ostrowski and Burdick (1998) that can predict motion from periodic shape-space changes for robotic systems could also be valuable.