In early work parameterizing snake kinematics, Shigeo (1994) proposed the serpenoid curve, wherein sinusoids approximate the shape of an undulating snake. The desired relative joint angle ϕ i ∗ ( t ) for joint i at time t follows Equation 1 with joint amplitude α h , frequency ω h , and phase offset between joints δ h . Biasing the relative joint angles with the steering term ϕ 0 induces whole body turning. ϕ i ∗ = α h s i n ω h t + δ h i − 1 + ϕ 0 (1)

Hirose’s model has become popular for snake robots due to its simplicity and versatility. The serpenoid curve has been used for both terrestrial and aquatic lateral undulation (Shigeo, 1994; Ostrowski and Burdick, 1996; Sato et al., 2002; Kelasidi et al., 2014). It has also been adapted for sidewinding (Tesch et al., 2009), concertina (Chen and Roth, 2023), and rectilinear (Tesch et al., 2009) gaits. Additionally, several authors have augmented the serpenoid curve to support obstacle aided locomotion (Rollinson and Choset, 2013; Travers et al., 2018; Travers et al., 2015) and mitigate slipping (Dehghani and Mahjoob, 2009). In our work, we use it as a basis for comparison with our SHC controllers. We select gait parameters α h = π 6 rads, frequency ω h = 5 π 9 rads/s, and δ h = 2 π 9 rads, consistent with Kelasidi et al. (2017).

To ensure locomotion along the x axis, we adapt the steering control in Kelasidi et al. (2017). The simulated robot heading ψ is computed as the average of all linkage angles with respect to the horizontal θ (See Figure 3; Equation 2). We assign a desired heading ψ ref using the steering constant Δ = 4 l (where 2 l is the length of one robot segment) and the center of mass displacement from the x axis, Y C M (Equation 3). Then, the steering term ϕ 0 is computed using proportional-integral (PI) control with integral gain k I , 0 = 0.2 π and proportional gain k p , 0 = 0.2 (Equation 4; Table 1). ψ = 1 n ∑ i = 1 n θ i (2) ψ ref = − a t a n Y C M Δ (3) ϕ 0 = − k I , 0 ∫ ψ ref − ψ d t + k p , 0 ψ ref − ψ (4)

The actuator torques u for each joint i are computed as in Kelasidi et al. (2014) with proportional gain k p = 50 , damping k d = 0.5 , and relative joint angles ϕ (Equation 5). u i = k p ϕ i ∗ − ϕ i + k d ϕ ̇ i ∗ − ϕ ̇ i (5)

Trajectory planning is enhanced by considering the interactions between the robot and environmental obstacles. Travers et al. (2015) define shape compliant control, where serpenoid curve parameters vary due to obstacle contact. Specifically, their robot senses externally applied torques on each joint and adjusts the reference joint amplitudes using an admittance controller, enabling the robot to brace against obstacles (Equation 6). M a ′ α h ̈ − α ̈ ref + D a ′ α h ̇ − α ̇ ref + K a ′ α h − α ref = τ ext ′ (6) In Travers et al. (2015), the admittance controller gains M a ′ , D a ′ , and K a ′ are functions of the robot configuration, and the reference amplitude α ref is held constant for a given joint. Their compliant control produces coordinated locomotion through an environment with regularly spaced pegs. Other methods of shape control have been explored for pipe climbing (Rollinson and Choset, 2013) and navigating rugged terrain (Travers et al., 2018).

We adopt a simplified shape controller based on Travers et al. (2015), where the robot narrows its gait by decreasing joint amplitudes during contact. Once contact ceases, the joint amplitudes return to the reference value α ref = π 6 (see Equation 7). This method reduces contact with the walls in confined spaces. We choose amplitude gain k α = 1500 , amplitude damping c α = 1050 , and tactile input gain k T = 500 . Tactile input z T = 1 if any linkages contact an obstacle and z T = 0 otherwise. Joints are uniformly impacted by tactile inputs from all segments. α ̈ h + c α α ̇ h + k α α h − α ref = − k T z T (7)

Replicating the anguilliform gait with SHC control provides an opportunity to assess the scalability of SHC movement primitive trajectories. In Rouse and Daltorio (2021), the authors demonstrate that weight parameters for SHC movement primitives visually represent and are proportional to the desired trajectory for a simple kinematic model. Inspired by their work, we examine how gaits can be adapted by scaling the weights. Specifically, we assess weight scaling as a means of switching from the terrestrial lateral undulation gait to anguilliform locomotion. Kelasidi et al. (2014) approximate anguilliform locomotion by linearly scaling the joint amplitudes in their serpenoid controller as in Equations 15, 16. g i , n = n − i n + 1 (15) ϕ i ∗ = α h g i , n s i n ω h t + δ h i − 1 + ϕ 0 (16) to implement an anguilliform gait with SHC control, we apply the scaling relationship in Equation 15 to the weight matrix (Equation 17). We expect that scaling the weights will proportionally change the amplitude of the trajectories for each joint. w i , j a = w i , j g i , n (17)

CPGs, both in biological systems and robots, are most useful when they can be adapted by sensory information. Here, we establish a mechanism for incorporating tactile signals into SHC movement primitives, with the goal of improving speed of the simulated snake robot as it undulates through confined spaces. We introduce the term k T z T into the SHC canonical state equation, where k T = 0.55 (Equation18). We assign z T = 1 if any linkage contacts an obstacle and z T = 0 otherwise, as in Section 2.2.2. During contact, kernels progress more quickly such that joints have less time to accelerate and joint amplitudes decline. The increased kernel progression rate also increases joint oscillation frequency, leading to higher thrust (Huang et al., 2019) that may enable the robot to overcome friction if contact persists. Since all tactile sensors affect every kernel, phase differences between joints remain constant and coordinated locomotion is maintained. In future work, we may consider alternative SHC formulations that allow tactile sensors to impact joints independently. τ d x = x 1 + k T z T α T − x ρ d t + η d t N 0,1 T (18)

Robots, like animals, benefit from optimizing the speed and energy efficiency of locomotion. In (Rouse and Daltorio, 2021), the authors demonstrate that weights in SHC movement primitives can be learned to replicate a known trajectory for a purely kinematic system. Here, we extend their work, improving the robot’s performance by adjusting weights with consideration to the system dynamics. Beginning with the SHC controller for lateral undulation (see Section 2.2.5), we simulate swimming for t max = 10 seconds in an obstacle-free environment. The forward velocity of the center of mass (Equation 19) and energy efficiency of locomotion are assessed. We define efficiency as the overall cost of transport (COT) (Equation 20), where t start = 0 s , G = 9.81 m s s is the acceleration due to gravity and m total is the total mass of the robot. v avg = 1 t max ∑ t = t start t max X ̇ C M d t t (19) C O T = ∑ t = t start t max ∑ i = 1 N joints u i , t ϕ ̇ i , t G m total X C M t max (20)

We establish a cost function (Equation 21) rewarding velocity, while penalizing poor efficiency. Note that the velocity is larger than the COT for lateral undulation such that this cost function prioritizes speed. c o s t = 5000 ‖ C O T ‖ − 5000 v avg (21) we numerically approximate the gradient of the cost function with respect to λ 1 − 4 , i for each joint i . The weights are improved via gradient descent with a learning rate η learn = 0.01 until the rate of change of weights d λ 1 − 4 , i ≤ 0.01 , requiring 155 iterations (Equation 22). The selected learning rate allows smooth improvement in performance, and changes become relatively small by the time the stop condition is reached (see Supplementary Figure). While our experiments in Sections 2.3.2, 2.3.3 include obstacles, the weights are not relearned for these conditions. The tactile sensing parameters are kept the same as in Section 2.2.5. w = w − η learn ⋅ ∇ c o s t w (22)

We first consider a simulated robot quiescent water with no obstacles. We adopt the framework proposed by Kelasidi et al. (2017), which addresses a submerged planar robot comprised of discrete linkages with elliptical cross sections. Fluid is viscous, incompressible, and irrotational in the reference frame of the robot. The authors develop equations of motion for the center of mass ( X C M , Y X M ) and linkage orientation θ including linear drag, nonlinear drag, and added mass in an unobstructed environment. Relative joint angles are derived from kinematics as ϕ i = θ i + 1 − θ i . We consider a robot with n = 6 segments, where each segment has a major radius e 1 = 0.09 m , minor radius e 2 = 0.035 m , half-length of l = 0.09   m , and a mass of m = 0.5   k g , except the head segment. The head segment is extended ( l 6 = 0.11   m and m 6 = 0.94   k g ) to accommodate electronic components for a physical robot, which is under development. Fluid parameters and fluid body interaction constants are kept consistent with Kelasidi et al. (2017) (fluid density ρ fluid = 1000   k g / m 3 , added mass coefficient C A = 1 , added inertia coefficient C M = 1 , drag coefficient along link axis C f = 0.03 , and drag coefficient perpendicular to link C D = 2). All simulation parameters can be found in Table 3.

The experiment is conducted for the serpenoid, lateral undulation SHC, anguilliform SHC, and learned SHC controllers. The simulation runs for 60 s, recording the robot’s position, joint torques, and joint velocities. Power consumption is approximated as P t = ∑ i = 1 n − 1 u i , t ϕ ̇ i , t . Then, we compare the swimming speed and COT for each controller. Additionally, we assess the differences in behavior for SHC controllers by examining the joint trajectories, kernel activity, and forcing term f ( x ) . At time t = 0 , the robot is at rest with all linkage angles θ i = 0 and the center of mass located at the origin.

Simulating the robot in structured environments allows us to assess our method for integrating tactile information with SHC control. We examine straight passages ranging from W = 0.2   m to W = 0.7   m in width, at intervals of 0.01   m . The minimum width allows 0.01   m of freedom on each side of the robot, while the maximum channel width permits the joints to oscillate freely if the robot is centered. The robot begins entirely in the channel with θ i = 0 and no contact with the walls. The overall COT and average velocity are calculated after 30 s of swimming for each case. This experiment addresses the serpenoid and lateral undulation SHC without sensing, as well as the serpenoid, lateral undulation SHC, and learned SHC with tactile sensing (+T).

Obstacles are treated as points connected by flat surfaces with stiffness k wall = 30,000 N m and damping c wall = 500 N m ⋅ s . We approximate each robot segment as a set of points spaced d x segment = 0.001   m apart and assume that points colliding with an obstacle penetrated the nearest surface. The normal distance d N and normal velocity v N between the point and surface are computed, then the normal force centered at the point is approximated as F N = ( k wall ⋅ d N + c wall ⋅ v N ) ⋅ d x segment . We incorporate friction via a linear model with stiction (Marques et al., 2016), where F S = μ s F N , F C = μ k F N , μ s = 0.8 , and μ k = 0.6 . The details are replicated in Equation 23. The variables v T , v 0 = 0.001 m s , and v 1 = 0.002 m s represent the tangential velocity, velocity at maximum static friction, and velocity at which only kinetic friction occurs (Table 3). Torque is obtained by integrating contact forces over the displacement from the link center of mass. The effects of obstacles on fluid behavior are ignored for simplicity, and fluid far from the robot remains quiescent. F f = v T v 0 F S sgn v T if  v T ≤ v 0 F S − v T − v 0 v 1 − v 0 F S − F C sgn v T if  v 0 < v T < v 1 F C sgn v T if  v T ≥ v 1 (23)

One of the most daunting aspects for exploratory robots is maintaining productive locomotion in complex terrain. We simulate the robot moving through channels with randomly generated feature (see examples in Figure 4). To establish consistent initial conditions, each channel begins with a 1.5 m long and 0.5 m wide straight section. The robot begins completely straight, at rest, and inside the channel. After the straight section, the upper surface of the channel is constructed such that x wall,i = x wall,i−1 + 0.025 and y upper,i = y upper,i−1 + N μ = 0 , σ = 0.03 . Points along the bottom surface y lower,i are computed likewise. If the distance between y upper,i or y lower,i and any point on the opposite boundary is less than the minimum channel width W min = 0.30 m ( 1.7 times robot width), the i t h upper and lower points are replaced. Similarly, the i t h points are regenerated if the distance between y upper,i and y lower,i exceeds W max = 0.5 m (2.8 times robot width) apart so that results do not simply reflect obstacle-free behavior. Due to random walk, the center lines shift laterally up to 1.6–3.0 robot widths over a single length of the robot, depending on the channel.

We simulate 60 s of swimming in each channel using the serpenoid and lateral undulation SHC controllers with and without tactile sensing. The learned SHC controller (+T) is also examined. Steering is disabled since the direction of the channel varies. We estimate the minimum number of channels required n channels,min = 61 using power analysis conducted in the software G*Power for 1-way ANOVA with effect size f p = 0.25 , power α p = 0.95 , and five groups. To account for invalid results, we increase the sample size to n channels = 65 . Specifically, data is excluded if the robot does not enter and remain within the uneven portion of the channel. We compute the COT and velocity using only time steps after the robot has fully entered the uneven section. Since the data does not follow normality assumptions, we assess statistical significance with Kruskal–Wallis testing, followed by Dunn-Šídák post hoc analysis. Furthermore we categorize the progress for each controller as continual, ceased, or intermittent. We specify that progress has ceased if the robot does not move forward by more than 0.1 m for at least 7.2 s (two gait cycles under unimpeded locomotion). We define intermittent progress where forward motion ceases for at least two gait cycles, but later resumes. The rates of continual, ceased, and intermittent progression for each controller are assessed.

In quiescent fluid with no obstacles, the SHC controller successfully replicates the lateral undulation gait. The SHC controller produces similar joint oscillations to the serpenoid controller, but with a smoother transition from rest to steady state oscillation (Figure 5A). For both methods, the heading naturally oscillates during undulation. The steering controllers compensate by increasing or decreasing the desired joint trajectories, causing joint amplitudes to differ from the reference value of α h = 30 ° .

Since the serpenoid and lateral undulation SHC controllers generate similar trajectories, they achieve comparable performance. The average velocity of the center of mass for both controllers is v avg,serpenoid ≈ v avg,LUSHC = 0.28 m s . The progression of the robot for all controllers is depicted in Figure 5 (also see Supplementary Files for Supplementary Video recording). (B) Due to differences in steering controller implementations, the SHC controller produces slightly better steady state efficiency ( C O T t start = 1 , L U S H C = 0.087 ) than the serpenoid controller ( C O T t start = 1 , s e r p e n o i d = 0.097 ) (Table 4).

A benefit of SHC control over serpenoid control is apparent when assessing the overall COT including transient effects ( t start = 0 ) . The conventional serpenoid curve yields a desired trajectory with nonzero initial joint angles and velocities, whereas most robots realistically start from rest. High accelerations, large joint torques, and high energy consumption are required to reduce the trajectory error, leading to an overall COT of 0.126 ( 30 % above steady state value). Shape-based control of the joint amplitude avoids this issue, albeit at increased complexity (Travers et al., 2015; 2018; Rollinson and Choset, 2013). In contrast, the SHC framework natively computes desired accelerations, then the desired velocity and position are evaluated by integration. This yields smooth kinematic trajectories such that the lateral undulation SHC controller achieves an overall COT only 1 % greater than its steady state value (Figure 5A).

By linearly scaling the SHC weights, we transform the lateral undulation gait to an anguilliform gait, with reduced anterior joint amplitudes. For example, we scale the weights for the anterior-most joint by g 5 = 14 % (Equation 15) relative to weights for the lateral undulation SHC controller. Consequently, the anguilliform SHC controller forcing term f 5 ( x ) magnitude is lower in comparison with the lateral undulation SHC (Figure 6A). This translates to a smaller joint amplitude (Figure 6C). The lateral undulation SHC controller exhibits a joint amplitude of 29.5 ° , while the anguilliform SHC controller’s anterior joint amplitude reduces to 4.0 ° ( 13.6 % ) , showing good correspondence with the scaling factor g 5 .

Scaling effects are more thoroughly characterized by considering the amplitudes of weights and ϕ for all joints (Figures 6B, D). The SHC weight amplitude for joint j is defined by the maximum value in the j t h row of the weight matrix (Figure 7). Anguilliform SHC weight amplitudes decrease linearly from the tail to the head, with similar trends emerging in joint amplitudes. The visible correspondence between weights, forcing term f ( x ) , and joint trajectories makes SHC MPs a transparent framework that robot operators can intuitively modify. Because lateral displacement is reduced, the anguilliform gait produces less thrust and slower swimming ( 0.166 m s ) than the lateral undulation gait. However, the anguilliform gait is more efficient by 26.1% ( C O T = 0.0718 ) relative to serpenoid control ( C O T = 0.0972 ) and 17% relative to lateral undulation SHC control ( C O T = 0.0867 ) . While we only simulate one sample and recommend more comprehensive experiments with hardware, results qualitatively agree with prior work showing that efficiency is optimized at lower amplitudes compared to velocity (Anastasiadis et al., 2023).

Our gait optimization for SHC control modifies the weight parameters to enhance speed of locomotion, while minimally affecting energy efficiency. The learned SHC controller achieves a 12.3 % higher overall velocity ( 0.32   m / s ) compared with the lateral undulation SHC controller. The steady state efficiency degrades, with the COT increasing from 0.087 in the lateral undulation gait to 0.10 in the learned gait.

The improved speed of locomotion in the learned SHC controller corresponds to changes in the controller weights and joint amplitudes. Relative to the lateral undulation SHC controller, the learned SHC controller possesses higher weight amplitudes for the posterior joints and lower weights for the anterior joints (7 (B)). The learned SHC weight amplitudes correspond loosely to the joint amplitudes (Figure 6D), which show a gradual decrease towards the head of the robot. Compared with anguilliform and lateral undulation SHC controllers, the learned controller shows weaker correspondence between weights and joint amplitudes because the distribution of weights, and thus acceleration patterns, may differ between joints.

When optimizing weight parameters for the learned SHC controller, the distribution of weights for each joint changes along with the peak amplitude. We compute the phases of weight distributions δ j w for each of the j joints using a centroid method (Equations 24–27). δ j , k = 360 k − 1 K 2 , 1 ≤ k ≤ K 2 (24) x j w = ∑ k = 1 K / 2 w j , k c o s δ j , k ∑ k = 1 K / 2 w j , k (25) y j w = ∑ k = 1 K / 2 w j , k s i n δ j , k ∑ k = 1 K / 2 w j , k (26) δ j w = a t a n y j w x j w (27)

Gait optimization produces weight phase offsets between joints Δ δ j w = δ j w − δ j − 1 w mostly near 30 ° , lower than the design value of δ h = 40 ° applied for other controllers (Figure 6E). The phase offset for the anterior joint pair may differ due to the longer segment length and increased mass, or because the Δ δ 5 w has little impact on the optimization cost function. For all controllers, the weight phase offsets identified in the weight matrix are good predictors for phase offsets in the joint trajectories Δ δ j = δ j − δ j − 1 (Figure 6F). Some deviation occurs due to steering and fluid effects in all controllers.