{'fraction_non_alphanumeric': 0.0349184961219929, 'fraction_numerical': 0.008150387800709872, 'mean_word_length': 4.387959274015051, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 23, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Reinforcement Learning of a CPG-regulated Locomotion Controller for a Soft Snake Robot

In this work, we present a learning-based goal-tracking control method for soft robot snakes. Inspired by biological snakes, our controller is composed of two key modules: A reinforcement learning (RL) module for learning goal-tracking behaviors given stochastic dynamics of the soft snake robot, and a central pattern generator (CPG) system with the Matsuoka oscillators for generating stable and diverse locomotion patterns. Based on the proposed framework, we comprehensively discuss the maneuverability of the soft snake robot, including steering and speed control during its serpentine locomotion. Such maneuverability can be mapped into the control of oscillation patterns of the CPG system. Through theoretical analysis of the oscillating properties of the Matsuoka CPG system, this work shows that the key to realizing the free mobility of our soft snake robot is to properly constrain and control certain coefficients of the Matsuoka CPG system, including the tonic inputs and the frequency ratio. Based on this analysis, we systematically formulate the controllable coefficients of the CPG system for the RL agent to operate. With experimental validation, we show that our control policy learned in the simulated environment can be directly applied to control our real snake robot to perform goal-tracking tasks, regardless of the physical environment gap between simulation and the real world. The experiment results also show that our method's adaptability and robustness to the sim-to-real transition are significantly improved compared to our previous approach and a baseline RL method (PPO).


I. INTRODUCTION
Due to their flexible geometric shapes and deformable materials, soft continuum robots have great potential to perform tasks in dangerous and cluttered environments, including natural disaster relief and pipe inspection [1]. However, planning and control of such robots are challenging, as these robots have infinitely many degrees of freedom in their body links, and soft actuators with stochastic, unknown dynamics and delayed responses.
We aim to develop a bio-inspired locomotion controller for soft robot snakes to achieve serpentine-like locomotion for setpoint tracking tasks. Specifically, we consider exploiting the properties of Central Pattern Generators (CPGs), which consists of a special group of neural circuits that are able to generate rhythmic and non-rhythmic activities for organ contractions and body movements in animals. Such activities can be activated, modulated, and reset by neuronal signals mainly from two directions: bottom-up ascendant feedback information from afferent sensory neurons or top-down descendant signals from high-level modules including mesencephalic locomotor region (MLR) [2] and motor cortex [3], [4].  In literature, bio-inspired control methods have been studied for the control design of rigid robots' locomotion, including bipedal [5]- [8] and serpentine locomotion [9]- [14]. The general approach is to generate motion patterns mimicking animals' behaviors and then track these trajectories with a closedloop control design. In [2], the authors developed a trajectory generator for a rigid salamander robot using Kuramoto CPGs and used low-level PD controllers to track the desired motion trajectories generated by the oscillator. Ryu et al. [11] established the velocity control CPG by adapting its frequency parameter with additional linear dynamics. In [13], the authors introduced a control loop that adjusts the oscillation patterns including frequency, amplitude, and phase of the oscillation for adapting to the changes in the terrain. Their result shows the advantage of the Hopf oscillator on the direct access to the oscillation patterns for different locomotion purposes. However, this work has not provided a way to maneuver the oscillation patterns autonomously. In [15], the Matsuoka oscillator is combined with the amplitude modulation method to realize steering control of a rigid snake robot. Sartoretti et al. [16] proposed a decentralized approach, where each actuator of an articulated rigid snake robot is controlled independently by a Neural Network (NN) controller learned with an end-toend reinforcement learning (RL) algorithm. Another recent work [14] employed a Spiking Neural Net (SNN) under the regulation of Reward-Modulated Spike-Timing-Dependent Plasticity (R-STDP) to map visual information into wave parameters of a phase-amplitude CPG net, which generates desired oscillating patterns to locomote a rigid snake robot chasing a red ball.
For a rigid snake robot, existing literature [2] has introduced a model-based control design combined with CPG for motion planning. We hypothesize that the proposed method should apply to a rigid snake robot as well. However, classical modelbased control for rigid snakes may not perform well for soft snake robots, due to the stochastic and difficult-to-model dynamics, as well as nonlinearity in the actuators. In previous work on a similar soft snake robot [17], the authors showed that the model-based approach results in around 20% error between real and estimation on a turning task (see ( Figure  10, [17])). In our control design, the bio-inspired approach is directly derived from the relation between CPG and rhythmic motion in the soft body of a real snake (see [18] for a detailed discussion of snake bio-mechanics.) Despite the success of bio-inspired control with rigid snake robots, the same control scheme may not work as desired for soft snake robots. This is because, in these approaches, the trajectories generated by CPG require high-performance lowlevel controllers for tracking. The tracking performance cannot be reproduced with soft snake robots due to the nonlinear, delayed, and stochastic dynamical response from the soft actuators.
To this end, we develop a bio-inspired learning-based control framework for soft snake robots with two key components: To achieve intelligent and robust goal-tracking with changing goals, we use model-free RL [19], [20] to map the feedback of soft actuators and the goal location, into control commands of a CPG network. The CPG network consists of coupled Matsuoka oscillators [21]. It acts as a low-level motion controller to generate actuation inputs directly to the soft snake robots for achieving smooth and diverse motion patterns. The two networks form a variant of cascade control with only one outer-loop, as illustrated in Fig. 1.
To reduce the time and samples required for learningbased control, we showed that several dynamic properties of the Matsuoka oscillators can be leveraged in designing the interconnection between two networks. We've proved that steering control can be realized by modulating both the amplitudes and duty cycles of the neural stimuli inputs of the CPG network, and velocity control can be realized by tuning the oscillating frequencies of the CPG net. These findings enable us to flexibly control the slithering locomotion with a CPG network given state feedback from the soft snake robot and the control objective.
This paper is an extension of our preliminary work [22] that designs a learning-based set-point tracking control for soft snake robots. In comparison to that, we make the following improvements: • Theoretical proof of steering maneuverability: We analyze the property of the biased oscillation in the Matsuoka oscillator. Using describing function analysis, we show that when the tonic inputs of the Matsuoka oscillator are bounded with certain constraints, the bias of the output signal becomes linearly related to the tonic inputs. This feature makes the steering control of the snake robot easier to learn for an RL agent. • Free-response Oscillation Constraints (FOC) for Simto-real: We investigate the transient property of the Matsuoka Oscillator from free-response oscillation to forcedresponse oscillation. Using this property, we introduce a fixed free-response tonic input signal to help regulate the amplitude and oscillation frequency of the forced tonic inputs which are generated by the RL policy. The new approach is referred to as Free-response Oscillation Constrained Proximal Policy Optimization Option-Critics with Central Pattern Generator (FOC-PPOC-CPG). This approach improves the transferability of the RL control policy learned in the simulation to the real robot.
• Improve reward density with potential field function: We newly introduce a potential-field-based reward shaping to accelerate the learning convergence. • Comprehensive Sim-to-real tests and analysis: We added new experiments comparing the learning efficiency and adaptability of the policy between the proposed method and vanilla PPO. Based on the experimental result for both simulation and reality, we show that our soft snake robot equipped with a properly designed "vertebrate" (the CPG system) can be more easily controlled by the RL agent. Our approach also achieves more reliable locomotion performance under various goal-reaching locomotion tasks that are unseen during the training process. The paper is structured as follows: Section II provides an overview of the robotic system and the state space representation. Section III presents the design and configuration of the CPG network. Section IV discusses the key properties of the CPG network and the relation to the design of an artificial neural network for the RL module. Section V introduces a curriculum and reward design for learning goal-tracking locomotion with a soft snake robot. Section VI presents the experimental validation and evaluation of the controller in both simulated and real snake robots. Our soft snake robot consists of n pneumatically actuated soft links. Each soft link of the robot is made of Ecoflex™ 00-30 silicone rubber. The links are connected through rigid bodies enclosing the electronic components that are necessary to control the snake robot. In addition, the rigid body components have a pair of one-direction wheels to model the anisotropic friction of real snakes. Only one chamber on each link is active (pressurized) at a time.

II. SYSTEM OVERVIEW OF THE SOFT SNAKE ROBOT
The configuration of the robot is shown in Figure 2. At time t, state h(t) ∈ R 2 is the planar Cartesian position of the snake head, ρ g (t) ∈ R is the distance from h(t) to the goal position, d g (t) ∈ R is the distance traveled along the goal-direction from the initial head position h(0), v(t) ∈ R is the instantaneous planar velocity of the robot, and v g (t) ∈ R is the projection of this velocity vector to the goal-direction, θ g (t) is the angular deviation between the goal-direction and the velocity direction of the snake robot. According to [23], the bending curvature of each body link at time t is computed by κ i (t) = δi(t) li(t) , for i = 1, . . . , 4, where δ i (t) and l i (t) are the relative bending angle and the length of the middle line of the i-th soft body link.
In [24], we developed a physics-based simulator that models the inflation and deflation of the air chamber and the resulting deformation of the soft bodies with tetrahedral finite elements. The simulator runs in real time using GPU. We use the simulator for learning the locomotion controller in the soft snake robot, and then apply the learned controller to the real robot.

III. DESIGN OF A CPG NETWORK FOR THE SOFT SNAKE ROBOT LOCOMOTION
In this section, we introduce our CPG network design consisting of interconnected Matsuoka oscillators [25], [26]. Primitive Matsuoka CPG: A primitive Matsuoka CPG consists of a pair of mutually inhibited neuron models. The dynamical model of the primitive Matsuoka CPG is given as follows: w ji y e j + u e i + c, where the subscripts e and f represent variables related to the extensor neuron and flexor neuron, respectively. The tuple (x q i , y q i ), q ∈ {e, f } represents the activation state and selfinhibitory state of i-th neuron respectively, z q i = g(x q i ) = max(0, x q i ) is the output of i-th neuron, b ∈ R is a weight parameter, u e i , u f i are the forced tonic inputs to the oscillator, and K f ∈ R is the frequency ratio. The set of parameters in the system includes the discharge rate τ r ∈ R, the adaptation rate τ a ∈ R, the mutual inhibition weights between flexor and extensor a ∈ R, the inhibition weight w ji ∈ R representing the coupling strength with the neighboring primitive oscillator, and the free-response oscillation tonic input c ∈ R (c = 0 in our previous work [22]). In our system, all coupled signals including x q i , y q i and z q i (q ∈ {e, f }) are inhibiting signals (negatively weighted), and only the tonic inputs are activating signals (positively weighted). In the current system, we have N = 4 primitive Matsuoka CPGs. For simplicity, we introduce a vector to represent all tonic inputs to the CPG net. Structure of the Matsuoka CPG Network for the Soft Snake Robot: Extending from a primitive Matsuoka CPG system to the multi-linked snake robot, we construct a CPG network shown on the right of Fig. 3. The network includes four linearly coupled primitive Matsuoka oscillators. It is an inverted, double-sized version of Network VIII introduced in Matsuoka's paper [25]. The network includes four pairs of nodes. Each pair of nodes (e.g., the two nodes colored green/yellow) in a row represents a primitive Matsuoka CPG (1). The edges correspond to the coupling relations among the nodes. In this graph, all the edges with hollowed endpoints are positive activating signals, while the others with solid endpoints are negative inhibiting signals. The oscillators are numbered 1 to 4 from head to tail of the robot. In order to build the connection between the CPG network and robot actuators, we define the output of the i-th primitive Matsuoka CPG as where A z represents the amplifying ratio of the difference between z e i and z f i of the i-th primitive oscillator. Given the Bounded Input Bounded Output (BIBO) stability of the Matsuoka CPG net [21], the outputs ψ = [ψ 1 , ψ 2 , ψ 3 , ψ 4 ] T from the primitive oscillators can be limited within the real region [−1, 1] by adjusting the ratio A z . We let ψ i = 1 for the full inflation of the i-th extensor actuator and zero inflation of the i-th flexor actuator, and vice versa for ψ i = −1. The actual pressure input to the i-th chamber is λ i · ψ i , where λ i is the maximal pressure input of each actuator. The primitive oscillator with green nodes controls the oscillation of the head joint. This head oscillator also contributes as a rhythm initiator in the oscillating system, followed by the rest parts oscillating with different phase delays in sequence. Figure 3 shows all activating signals to the CPG network. Configuring the Matsuoka CPG Network: To determine the hyper-parameters in the CPG network that generate a more efficient locomotion pattern, we employed a Genetic Programming (GP) algorithm similar to [27]. In this step, all tonic inputs are set to be 1 for the simplicity of fitness evaluation.
We define the fitness function-the optimization criteria-in GP as F (v g,T , θ g,T , d g,T ) = a 1 |v g,T | − a 2 |θ g,T | + a 3 |d g,T |, where g indicates a fixed goal initiated in the heading direction of the snake robot, T indicates the terminating time of fitness evaluation for each trial, and all coefficients a 1 , a 2 , a 3 ∈ R + are constants 1 .
To achieve stable and synchronized oscillations of the whole system, the following constraint must be satisfied [21]: where τ a , τ r , b > 0. To satisfy this constraint, we can set the value of b much greater than both τ r and τ a , or make the absolute difference |τ r − τ a | sufficiently small. In other wordss, this fitness function is a weighted sum over the robot's instantaneous velocity, angular deviation, and total traveled distance on a fixed straight line at terminating time T . Under this scenario, a better-fitted configuration is supposed to maintain oscillating locomotion and reaches faster locomotion speed |v g,T | along the original heading direction at termination T . In addition, the locomotion pattern is also required to have a less angular deviation from the initial heading direction (with a small |θ g |), and with overall a longer traveled distance along the initial heading direction (|d g |).
The desired parameter configuration found by GP is given in Table. II in the Appendix.

IV. MANEUVERABILITY ANALYSIS AND DESIGN OF THE LEARNING-BASED CONTROLLER WITH THE MATSUOKA CPG NETWORK
When provided with equally constant tonic inputs, the designed Matsuoka CPG system can generate stable oscillation patterns to efficiently drive the soft snake robot slithering forward. However, a single CPG network cannot achieve intelligent locomotion and goal-tracking behaviors with potentially time-varying goals by itself. Free turning and accelerating (or decelerating) are the fundamental skills for an artificial agent to realize autonomous locomotion for goal-tracking tasks. In this paper, we denote these two maneuverability demands as -steering control and speed control (see Fig. 4). The later parts will focus on investigating the properties of the Matsuoka CPG system to prove that it is controllable from both steering and speed control perspectives. We will design a proper connection between RL actions and controllable coefficients of the Matsuoka CPG system so that both steering and speed control of the snake robot can be efficiently learned by the RL agent.
For steering control, we will prove that the bias of tonic inputs is linearly proportional to the bias of the CPG output in both amplitude and duty cycle dimensions. This property inspires a rule that transforms the action outputs of the RL policy into the tonic inputs of the CPG system.
Next, we excavate two mechanisms that are helpful for speed control. First, we show that the frequency ratio coefficient K f allows the RL agent to tune the locomotion speed by directly adjusting the oscillation frequency. Secondly, by introducing the free-response oscillation constraint, we provide a way to adjust the converging amplitude of the oscillation driven by the RL agent. With experiments, we will show that the freeresponse oscillation constraint is very helpful for reducing performance drop in the sim-to-real problem.

A. Steering control with imbalanced tonic inputs
Most existing methods based on CPG realize steering by either directly adding a displacement [2] to the output of the CPG system, or using a secondary system such as an artificial neural network to compose the weighted outputs from multiple CPG systems [5]. In this section, we present a different approach based on the maneuverability of the Matsuoka oscillator-tuning tonic inputs to realize the biased wave patterns of CPG outputs for steering the slithering locomotion of the soft snake robot 2 .
For the RL controller to steer our snake robot smoothly through the Matsuoka CPG system, we need to find a clever way to make the steering dynamics easy to learn for the RL algorithm. In other words, the relation between tonic inputs and the output bias of each primitive Matsuoka oscillator in the CPG network needs to be simple and clear. The winning secret hides in the relationship between each pair of tonic inputs u e and u f in a primitive Matsuoka oscillator. In the original design of the Matsuoka oscillator, the flexor and extensor tonic inputs are independent of each other. This setting not only increase the dimension of action space for the RL agent but also makes the relationship between tonic inputs and the output bias more complicated. To simplify this problem, we first introduce a new relation named exclusiveness to reform the relation between u e and u f . Definition 1. (Exclusiveness) For two real signals u(t) and v(t), and a known bounded range D : [a, b] where D ⊆ R, we say u(t) and v(t) are exclusive to each other in range D when From our previous work [22], we have statistically observed that the steering bias of a primitive Matsuoka oscillator is proportional to the amplitude of u e when u e and u f are exclusive within the range [0, 1]. This key observation allows us to apply a dimension reduction technique to the input space of the CPG net: Instead of controlling u e i , u f i for i = 1, . . . , n for a n-link snake robot, we only need to control u e i for i = 1, . . . , n and let u f i = 1 − u e i . As the tonic inputs have to be positive in Matsuoka oscillators, we define a four dimensional action vector a = [a 1 , a 2 , a 3 , a 4 ] T ∈ R 4 and map a to tonic input vector u as follows, This mapping bounds the tonic input within [0, 1]. The reduced input dimension enables a more efficient policy search in RL. Based on this design, we show that there are certain combinations of tonic inputs in a Matsuoka oscillator that are capable of generating imbalanced output trajectories and therefore result in the turning behavior of the robot. We present three possible cases of the forced tonic inputs that could maneuver the turning behavior of the snake robot: • The two tonic inputs are different constants. • The two tonic inputs are wave functions with imbalanced duty cycles. • The two tonic inputs are wave functions with imbalanced duty cycles, and both wave functions are added by different constant offsets. It is noted that the third case is a linear combination of the first two. As a result, as long as the first two cases are proved to share the same property, the third one naturally holds. Next, we provide the frequency domain analysis of the Matsuoka oscillator to explain why the first two cases of tonic inputs share the same property that enables imbalanced oscillation for the turning behavior. Steering with biased amplitude of constant tonic inputs: To show that a pair of tonic inputs with different amplitudes can result in a biased oscillating output trajectory, we need to find out the relation between the bias of the output z i and the tonic inputs when the tonic inputs are exclusive to each other. When studying this relation, we focus on the case where the oscillation converges to a stable limit cycle, which means that the system is harmonic. According to [26], given a harmonic primitive Matsuoka oscillator in Eq. (1), K n is the amplitude ratio between state x q i and output z q i for q ∈ {e, f }, which has the following form where τ r and τ a are the discharge rate and the adaptation rate in Eq. (1), and parameter a is the mutual inhibition weight between flexor and extensor of a primitive Matsuoka oscillator. The derivation of K n can be found in Appendix A-A. When both u e and u f are constants and are exclusive to each other in [0, 1], the forced response of the Matsuoka oscillator can be considered as zero and the system's output only consists of free-response oscillation. Under this condition, we have the following relation between the output oscillation bias bias(ψ) and the tonic input bias bias(u) (u = u e − u f ). Proposition 1. If a primitive Matsuoka oscillator has constant tonic inputs u e and u f being exclusive to each other, then the oscillation bias of its output signal ψ(t) satisfies where K n = (τ r + τ a )/(τ a a) is a constant coefficient.
In Proposition 1, Eq. (7) suggests that there is a linear relationship between the tonic input u q (q ∈ {e, f }) and the bias of the corresponding output given the constraint by Eq. (5) is satisfied by u e , u f . In addition to the proof, we validate this claim in the simulation of a single primitive Matsuoka oscillator. We calculate the mean oscillation bias (numerical average) of the simulated state output ψ and compare it with the estimated bias based on Eq. (7) (linear reference). Figure 5 shows the value of bias(ψ) given the input u e changing from 0 to 0.5 and u f = 1 − u e in a primitive Matsuoka oscillator. Here, we approximate the bias of the steady-state output oscillation trajectory by taking the time average of the trajectory 3 . From Fig. 5, we can observe a linear mapping between tonic input u e of a primitive oscillator and the bias of output when u e and u f are exclusive to each other. This result matches the prediction presented by Eq. (7).
In Fig. 5, the dashed line shows a threshold of bifurcation that separates the converging property of the Matsuoka oscillator into stable equilibrium and limit cycle oscillation. The system evolves in a periodical way only when the value of u e is above the threshold indicated by the dashed line. It is noticed that in all sub-figures of Fig. 5, the numerical average of ψ shows a bifurcation from fixed-point equilibrium to limit cycle equilibrium as the value of tonic input u e increases. The curve of the numerical average only matches the linear reference in the limit cycle region. This is because according to the perfect entrainment assumption [26], Proposition 1 only holds when the system is contracting to a harmonic limit cycle. Otherwise, when the system is converging to a set point, approximating its state trajectory with a sine wave is inaccurate. However, 3 Based on Fourier series analysis, given a continuous real-valued Pperiodic function ψ(t), the constant term of its Fourier series has the form 1 P P ψ(t)dt. a general model-free reinforcement learning algorithm often converges to a discontinuous control policy unless enforced by extra regularization techniques [28]. As a result, we need the features in Proposition 1 to hold even when the tonic input signals are oscillatory and discontinuous. In the next paragraph, we show that the biased duty cycle of the tonic input signal also has a linear relationship with the output bias of the Matsuoka oscillator.
Steering with the biased duty cycle of periodic tonic inputs: We consider a different approach to control the steering of the snake robot given that both u e i and u f i are pulse wave functions and are exclusive to each other. We define the duty cycle of a wave function as D(·). Given a Matsuoka oscillator in Equation (1), we aim to find out the relationship between the bias of ψ i and the duty cycle of the tonic inputs. Before going further, we first define a relation between two periodic signals named "entrainment" based on the related theory in [26], [29].
Definition 2. (Entrainment) Given a CPG system with its natural frequency ω n > 0. If a CPG system's output is synchronized to the forced input frequency ω, then this system is entrained with the forced input signal. The relation between the CPG output and the forced input signal is called entrainment.
According to the perfect entrainment assumption [26], let u i be resonant to ψ i . Let the period of ψ i be T = 2π (it is easy to prove that a different value of T would not change the following equation), based on the Fourier expansion, the bias of u i can be expressed as bias(u i ) = 1 2π Because the bias terms of x i and u i are time-invariant, from eq. (A.6), we can extract the bias component to form a new equation as follows bias(y i ) = bias(ψ i ), If a primitive Matsuoka oscillator has pulse wave tonic inputs u e and u f , and both tonic inputs are exclusive to each other, then the oscillation bias of the system output ψ(t) approximately satisfies where d(·) is the duty cycle function, K m and M are a constant coefficients.
(For details of the proof, please refer to Appendix B-B.) The simulated result also supports Proposition 2. Fig. 6 shows a linear relationship between the duty cycle and the output bias sampled from the state trajectory of a primitive Matsuoka oscillator. It is noticed that similar to Proposition 1, the second proposition only holds when the output state ψ converges to a stable limit cycle.
Since we have proved that both Proposition 1 and Proposition 2 are true based on Definition 1, we can conclude the following remark, Remark 1. If a primitive Matsuoka oscillator has periodical tonic input signals u e and u f that are exclusive to each other, then the oscillation bias of the system output is linearly related to the bias between u e and u f . Generally, the linear velocity of serpentine locomotion is affected by the snake's oscillation amplitude and frequency. In this subsection, we show that the amplitude and frequency can be controlled by two coefficients of the Matsuoka CPG system to change the locomotion speed of the soft snake robot.

B. Velocity control with frequency modulation
First, the following relation between the frequency ratio K f and the natural frequencyω i of the i-th oscillator is established in [26,Eq. (5), Eq. (6)], Second, the oscillating amplitudeÂ i of the i-th oscillator is linearly proportional to the amplitude of free-response oscillation tonic input c when c > 0 and u e i , u f i are constants [26], that is,Â where A(·) is the amplitude function of a rhythmic signal. Equations (11) and (12) show that the frequency and amplitude of the Matsuoka CPG system are independently influenced by the frequency ratio K f and the free-response oscillation tonic inputs c. Therefore, these two coefficients can be considered major factors for the Matsuoka CPG system to control the velocity of the soft snake robot's locomotion. In Fig. 7, we collect 2500 uniform samples within the region c ∈ [0.4, 0.8], and K f ∈ [0.45, 1.05] and record the velocities generated in the simulator. From Fig. 7, we observe that with a fixed c, the average velocity increase nearly monotonically with the frequency ratio K f . We also observe that with the same K f , the change of c does not affect the locomotion velocity significantly. But with a different value of c, the efficiency of K f in affecting the locomotion velocity is different. This means that we can mainly use K f to adjust the locomotion velocity, but the value of c needs to be carefully selected 4 . Given this analysis, we use K f to control the velocity of the robot. It is noted that the frequency ratio K f only influences the strength but not the direction of the vector field of the Matsuoka CPG system. Thus, modulating K f will not affect the stability of the whole CPG system.

C. Modulating forced-response oscillation amplitude with freeresponse oscillation tonic input constraint
In Eq. (12), we show that the free-response oscillation tonic input c could affect the output amplitude of the Matsuoka oscillator when u e and u f are constants. We will further discuss how a positive c could maintain the amplitude of u e and u f when they are periodic signals. In the experiment section, we will show that this feature of c can significantly improve the sim-to-real performance of our control framework.
In our previous work [22], we set c = 0. In that situation, there is no free-response oscillation in the system. When a Matsuoka oscillator has no free-response oscillation pattern, its output oscillation amplitude and bias are only controlled by the forced input signal given by the control tonic inputs u e and u f . When the inertia in the simulated learning environment is high and the contact friction force is low, the RL agent will learn to generate the forced-response oscillation tonic inputs with very small amplitude to keep a more stable heading direction during the locomotion. However, if we need the RL control policy to be able to initiate CPG oscillation with a larger amplitude on the real robot, the learned policy will not meet the requirement.
Inspired by Matsuoka's draft [29], the amplitude A of u e and u f (they share the same amplitude according to Eq. (5)) must satisfy the inequality A > A 0 (ω) to completely entrain with the output ψ of the Matsuoka oscillator (In [29], A 0 is the entrainment threshold for the amplitude of u e and u f to synchronize the output of the Matsuoka oscillator). Otherwise, if A ≤ A 0 (ω), the output pattern will tend to become a mixture of free-response oscillation and forced-response oscillation, which intervenes the turning behavior of the control policy during the goal-reaching tasks. The equation of A 0 (ω) is given as follows where A n > 0 is the free-response oscillation amplitude and ω n = 1 τrK f (τr+τa)b τaa − 1 is the free-response oscillation frequency [26]. The detailed calculation of A 0 is provided in Appendix. A-C. From the above equation, when c = 0, A 0 (ω) ≡ 0, there is no limiting threshold for the control policy to entrain the CPG output with u e , u f , regardless of the amplitude A of u e and u f . When c is assigned with a positive value, both A n (determined by c and u e i , u f i ) and ω n (determined by K f ) will make A 0 (ω) > 0. Remember that only when A > A 0 (ω) could the forced tonic inputs u e , u f entrain the CPG output ψ, the control policy will have to increase the amplitude of u e , u f to control the CPG system effectively. It is also noted that A 0 (ω) → 0 as ω → ω n , therefore there are two ways for the RL agent to realize the entrainment status: one is keeping the oscillation frequency ω close to the freeresponse oscillation frequency ω n , and the other is increasing the value of A to make A > A 0 . The combination of the two directions encourages the forced tonic input u e and u f to approach desired oscillating amplitude and frequency. Based on this special property of the Matsuoka oscillator, we propose a new method -FOC-PPOC-CPG to enforce better entrainment between RL control signals and the CPG states.
According to the property of A 0 (ω) and c, we can choose different values of c to keep the oscillation amplitude of the Matsuoka oscillator at different levels. One previous work [30] has shown that the oscillation amplitude of the Matsuoka oscillator can be used to adjust the serpentine locomotion velocity of a rigid snake robot in a different environment with different friction coefficients. This result inspires us to use c to adapt the locomotion speed of the soft snake robot to a different environment with various contact surfaces. If this holds true for most contact surfaces in practice, we will have a tool for the sim-to-real problem of our soft snake robot from the biomimic point of view, instead of fully relying on the environment-based approaches such as domain randomization or other data augmentation methods, which are computationally expensive.
In the later part of this paper, our experiment results (see Section VI-D) will show that this approach can naturally improve the sim-to-real performance by constraining the oscillating patterns of the Matsuoka CPG system.

D. The Neural Network Controller
We have now determined the encoded input vector of the CPG net to be vector a (tonic inputs) and frequency ratio K f . This input vector of the CPG is the output vector of the NN controller. The input to the NN controller is the state feedback of the robot, given by s = [ρ g ,ρ g , θ g ,θ g , κ 1 , κ 2 , κ 3 , κ 4 ] T ∈ R 8 (see Fig. 2). Next, we present the design of the NN controller.
The key insight for the design of the NN controller is that the robot needs not to change velocity very often for smooth locomotion. This means the updates for tonic inputs and the frequency ratio can be set to be at two different time scales.
With this insight, we adopt a hierarchical reinforcement learning method called the option framework [31], [32] to learn the optimal controller. The controller uses the tonic inputs as low-level primitive actions and frequency ratio as high-level options of the CPG net. The low-level primitive actions are computed at every time step. The high-level option changes infrequently. Specifically, each option is defined by I, π y : S → {y} × dom(a), β y where I = S is a set of initial states, and π y is the intra-option policy. By letting I = S, we allow the frequency ratio to be changed at any state in the system. Variable y ∈ dom(K f ) is a value of frequency ratio, and β y : S → [0, 1] is the termination function such that β y (s) is the probability of changing from the current frequency ratio to another frequency ratio.
The options share the same NN for their intro-option policies and the same NN for termination functions. However, these NNs for intro-option policies take different frequency ratios. The set of parameters to be learned by policy search include parameters for intra-option policy function approximation, parameters for termination function approximation, and parameters for high-level policy function approximation (for determining the next option/frequency ratio). Proximal Policy Optimization Option-Critics (PPOC) in the openAI Baselines [33] is employed as the policy search in the RL module.
Let's now review the control architecture in Figure 3. We have a Multi-layer perceptron (MLP) neural network with two hidden layers to approximate the optimal control policy that controls the inputs of the CPG net in Eq. (1). The output layer of MLP is composed of action a (green nodes), option in terms of frequency ratio (pink node), and the terminating probability (blue node) for that option. The input of NN consists of a state vector (yellow nodes) and its output from the last time step. The purpose of this design is to let the actor network learn the unknown dynamics of the system by tracking the past actions in one or multiple steps [5], [34], [35]. Given the Bounded Input Bounded Output (BIBO) stability of the Matsuoka CPG net [21] and that of the soft snake robots, we ensure that the closed-loop robot system with the PPOC-CPG controller is BIBO stable. Combining with Eq. (5) which enforces a limited range for all tonic inputs, this control scheme is guaranteed to generate bounded inputs, which lead to bounded outputs in the system.

V. CURRICULUM AND REWARD DESIGN FOR EFFICIENT LEARNING-BASED CONTROL
Here we introduce the design of the curriculum and reward function for learning a goal-tracking controller given the proposed PPOC-CPG scheme.

A. Task curriculum
Curriculum teaching [36] is used to accelerate motor skills learning given complex goal-tracking tasks. We design the curriculum such that the agent starts with easy-to-reach goals at level 0. As the level increases, the agent learns to perform more complex goal-tracking tasks.
Specifically, the curriculum levels are designed as follows: At the task-level, i, the center of the goal is sampled from the 2D space based on the current location and head direction of the robot. For each sampled goal, we say the robot reaches the goal if it is r i distance away from the goal. The sampling distribution is uniform in the fan area determined by the range of angle θ i and distance bound [ρ l i , ρ u i ] in the polar coordinate given by the predefined curriculum.
As shown in Fig. 8, when the task-level increases, we have ; that is, the robot has to be closer to the goal in order to succeed and receive a terminal reward, the goal is sampled in a range further from the initial position of the robot. We select discrete sets of {r i }, {θ i }, [ρ l i , ρ u i ] and determine a curriculum table. A detailed example of the learning curriculum is given in Table III. We train the robot in simulation starting from level 0. The task-level will be increased to level i + 1 from level i if the controller reaches the desired success rate σ i , for example, σ i = 0.9 indicates at least 90 successful completions of goal-reaching tasks out of n = 100 trials at level i.

B. Reward design
The design of the reward function is to guide the robot to the set point goals. We consider building the artificial potential field [37] such that the robot is attracted by the goal g. We use a simple conical potential field for each goal. For any position represented by coordinate x in Cartesian space, let vector e g = x g − x, the norm ||e g || indicates the distance between the position of the robot's head and the goal. The constant attracting force at x becomes f g (x) = e g ||e g || .
Considering the single goal-tracking scenario without obstacles, we have the potential field reward for goal-tracking as where v s is the velocity vector of the soft snake robot.
Combining with the definition of goal-reaching tasks and their corresponding level setups, the reward at every time step is defined as where c v , c g ∈ R + are constant weights, v g is the velocity towards the goal, θ g is the angular deviation from the center of the goal, r k defines the goal range in task-level k, for k = 0, . . . , i, ρ g is the linear distance between the head of the robot and the goal, and I(ρ g < r k ) is an indicator function that outputs one if the robot head is within the goal range for task-level k. This reward trades off two objectives. The first term, weighted by c v , encourages high locomotion velocity toward the goal. The second term, weighted by c g , rewards the learner based on the position of the robot to the goal, and the level of the curriculum the learner has achieved for the goal-reaching task. For every task, if the robot hasn't entered the goal range, it will receive a potential field reward only. When the robot enters the goal range in task-level i, it will receive a summation of rewards 1/r k for all k ≤ i (the closer to the goal the higher this summation), shaped by the approaching angle θ g (the closer the angle to zero, the higher the reward).
If the agent reaches the goal defined by the current task-level, a new goal is randomly sampled in the current or next level (if the current level is completed). There are two failing situations, where the desired goal will be re-sampled and updated. The first situation is starving, which happens when the robot stops moving for a certain amount of time. The second case is missing the goal, which happens when the robot keeps heading in the wrong direction as opposed to moving towards the goal (v g (t) being negative) for a certain amount of time. In this section, we evaluate the proposed method in both simulation and real environments. We first introduce the experimental setup to explain how the dynamic data of the robot is collected during locomotion, as well as the training configuration for the RL algorithm. Then we compare the properties of control signals between our method and the vanilla PPO as locomotion controllers for goal-reaching tasks in simulation. In the comparison analysis, we highlight the performance drop of each method from simulation to real to extrapolate the advantage of our method. Last, we further test and analyze the sim-to-real robustness of the up-to-date method, including sim-to-real performance in difficult goalreaching trajectories that are never seen at the training stage, and the real robot performance against disturbance during a goal-reaching process. All the experimental results provide promising evidence of our method as an effective and robust controller for intelligent soft snake locomotion.

A. Experimental Setup
Environment Sensing and Data Collection: The physical states of the real snake robot are captured by a single web  camera hanging on the ceiling of the experiment room. The robot body detection is realized by using Aruco -a library in OpenCV for QR codes detection and localization [38]. These QR codes are printed and attached to the rigid bodies of the snake robot and the goal position. Fig. 9 shows the experimental setup for the real snake goal-reaching tasks. Once the QR codes on the robot bodies and the goal(s) are detected, their pixel-wise coordinate vectors are calculated with distortion corrected. Given the camera calibration information, we can translate the pixel data of all QR codes into the real world 2D coordinates, and then transform it into positional information and the body posture of the robot. The control policy function running on the desktop computer receives the observation states and generate the control commands and passes them through WiFi communication. The ESP32 chips on the snake bodies translate the commands into Pulse Width Modulation (PWM) signals to activate or deactivate the valves [24], [39] on the snake robot. Reinforcement Learning Configuration: We use a fourlayered NN with 128 × 128 hidden layer neurons as a general configuration for the actor and critic networks of all RL methods mentioned in this section. The back-propagation of the critic net was done with Adam Optimizer and a step size of 5 × 10 −4 . For data collection of each trial trajectory, the starvation time for the failing condition is 60 ms. The missing goal criterion is triggered whenever v g (t) (the velocity on the goal-direction, see Fig. 2) stays negative for over 60 time steps. In order to compensate for mismatch in the simulation and real environment, most notably friction coefficients, we employ a domain randomization technique [40], in which a subset of physical parameters are sampled from several uniform distributions. The range of distributions of Domain Randomization (DR) parameters used for training are in Table  IV (see Appendix C).
Especially, for PPOC-CPG, we first train the policy net with fixed options (at this moment, the termination probability is always 0, and a fixed frequency ratio K f = 1.0 is used). When both the task-level and the reward cannot increase anymore, we allow the learning algorithm to change the option, i.e., pick a different frequency ratio K f along with termination function β, and keep training the policy until the highest level in the curriculum is passed.
A newly included parameter in this paper is the freeresponse oscillation tonic input c. In the PPOC-CPG method, the value of c is equivalently considered zero since it is not formally introduced in the previous paper [22]. According to the definition of A 0 (ω) (Eq. (13)), the amplitude of both u e and u f needs to be greater than A 0 (ω) in order to dominate in controlling the output pattern of the Matsuoka CPG system. The value of A 0 (ω) should not be greater than the upper bound of u e and u f , which is 1 defined by Eq. (5). Among a group of candidates ranging from 0.25 to 2, we choose c = 0.75 as our free-response oscillation constraint for the FOC-PPOC-CPG controller. This value is valid for our system because when we plug in c = 0.75 and all other coefficients of the CPG network (Table II) to Eq. (13), the result shows It is noted that the range of ω here is calculated from multiple sampled sequences of u e and u f recorded in the real snake reaching tasks. Since we are testing the sim-to-real performance, all methods involved in this comparison are trained in the simulator for sufficiently long iterations (12500 episodes) to ensure convergence. Each method is trained 10 times with different random seeds and the controller with the best performance is selected to be tested on the real robot. All curriculum parameters (Table III) and domain randomization parameters (Table IV) are fixed for all three methods involved.
The whole training process of each method runs on 4 simulated soft snake robots in parallel on a workstation equipped with an Intel Core i7-9700K, 32GB of RAM, and one NVIDIA RTX2080 Super GPU.

B. Verification of steering property of PPOC-CPG
We use a simulated experiment to show that our FOC-PPOC-CPG control policy has learned the turning behavior with the biased tonic input signal, and the Matsuoka CPG system can linearly map the biased tonic input to the biased actuation signal as Proposition 1 and Proposition 2 predicted. In the experiment, we test the converged FOC-PPOC-CPG policy on multiple set-point goals placed in certain directions (−90 • , −70 • , −60 • , −45 • , −30 • ) with a fixed distance (1 meter), which approximately represent the desired turning angles of the locomotion tasks. For each goal position, we run 5 trials and record the tonic inputs data and CPG output data of the head CPG node of the soft snake robot. The reason for choosing the head node is because this node's behavior best reflects the desired steering direction of the RL agent. Figure 11a shows a violin plot of the tonic input bias and the CPG output bias for different turning angles (the bias signals are calculated by Eq. (8)). It is observed from Fig. 11a that both bias signals are monotonically related to the desired turning angle (initial goaldirection). Figure 11b shows the linear regression result based on all data points. We can observe a clear linear relationship between bias(ψ) and bias(u) of the head CPG node (with the coefficient of determination equal to 0.939, a value closer to 1 indicate higher linearity). This result provides stronger support for Proposition 1 and Proposition 2.

C. Control signal comparison between PPOC-CPG and vanilla PPO
First, we compare PPOC-CPG and vanilla PPO in terms of the smoothness of the control input learned in simulation. We train both PPOC-CPG and vanilla PPO in the same environment until convergence. Figure 10 shows segments of the control signal ψ 1 generated by the vanilla PPO controller and PPOC-CPG controller controlling the simulated soft snake robot in a straight line goal-tracking task. From Fig. 10a it is observed that the signal generated by the vanilla PPO policy oscillates at a relatively higher frequency (about 10Hz on average) with irregular oscillation patterns. Such kind of control signals are not feasible for the actuators in reality. This is because the inflation and deflation of soft air chambers on the snake robot have a certain delay so that the soft pneumatic actuators are not able to track fast oscillating signals. On the other end, the curve in Fig. 10d shows that the agent trained with PPOC-CPG can converge to a stable limit cycle trajectory at a relatively lower but more natural frequency (1.6Hz) for serpentine locomotion. Our approach shows its advantage of being able to generate smoother oscillatory control signals even when the inputs to the CPG system are discontinuous. In the next section (Section VI-D), these methods are compared in the real robot to demonstrate the advantage of the proposed PPOC-CPG control. Fig. 10 also compares the phase plane portraits recorded at different time stages of the two learning methods. From Fig. 10b and Fig. 10c, we observe that the oscillating signal generated by vanilla PPO policy performs irregular oscillation in the first 300 time steps, and cannot converge to a stable limit cycle when it evolves to time step 700. While in Fig. 10e and Fig. 10f, despite a little deviation from the first 300 time steps, the outputs of the CPG network eventually converge to a stable limit cycle within 700 time steps. This experiment shows that the CPG system is capable of stabilizing the oscillation pattern in simple locomotion tasks for the soft snake robot.
D. Sim-to-real Performance of FOC-PPOC-CPG 1) Performance comparison with original PPOC-CPG and Vanilla PPO: Since FOC-PPOC-CPG is designed for improving the sim-to-real performance of the PPOC-CPG method, we first compare the sim-to-real performance of the FOC-PPOC-CPG with the original PPOC-CPG and Vanilla PPO in single goalreaching tasks. For the real robot tests, all three controllers trained by the simulator are directly applied without further training. To ensure convergence, we test the controllers by setting goals in three directions (mid, left and right) with fixed angles, distances, and an accuracy radius of r = 0.175 meters. Each direction takes 10 trials for all three methods in both simulation and reality.
To evaluate the sim-to-real performance, we calculate the average speed and the success rate for goal-reaching tasks collected from both simulation and real experiments. According to Section IV-C, the contact resistance forces in the simulator are smaller than in the real environment, when applying the RL control policy learned in the simulator directly to the real robot, the performance of the real robot is often worse than the simulated agent. In the rows of real robot evaluations in Table I, we use down-arrows and percentage values to show the extent of performance drop compared to the simulating performance with the same method. The evaluation results of each method take an entire column in the table. From Table I, although the Vanilla PPO controller learns the best locomotion speed in the simulator at the cost of goal-reaching accuracy, its locomotion pattern cannot fit the real robot well. It is observed that the real robot experiences a drastic drop in performance on both locomotion velocity (80.7%) and success rate (42.1%). For the original PPOC-CPG, though it has achieved an overall better performance than Vanilla PPO, its sim-to-real performance drop is still relatively high, with a 54% of speed drop and 17% of accuracy drop. After adding the free-response oscillation constraint to the CPG system, the new policy reaches almost the same performance as the   original PPOC-CPG in the simulator. In Section IV-C we have shown that the free-response oscillation tonic input c > 0 could help maintain the oscillation amplitude of the control signal of FOC-PPOC-CPG during the learning process. It is noticed that the maintained amplitude of the control signals does not improve the locomotion speed and goal-reaching accuracy at the training stage in the simulation. But when the learned policy of FOC-PPOC-CPG is applied to the real robot without further training, it performs significantly better than the previous two methods in both velocity and success rate. 2) Performance in reaching unseen goals: We also investigate the sim-to-real performance of FOC-PPOC-CPG in harder goal-reaching tasks. Figure 13 compares the head trajectories in Cartesian space for two different setups of way-point goals. The testing trajectories include a square turning trajectory for testing consecutive sharp turning in the same direction (Fig. 13b), and a zigzag trajectory for testing continuous sharp turning in opposite directions (Fig. 13a). Both way-point goal series have sharper turning angles than the highest level in the training curriculum in Table III. Video "Half square trajectory sim2real.mp4" and "Zigzag trajectory sim2real.mp4" provide the dynamic view of Fig. 13a and Fig. 13b respectively. From the example videos, it is observed that in both trajectories, the velocity drop of the real robot is still around 10%, which is not getting worse than single goal-reaching tasks in Table I. It is noted that in both Fig. 13a and Fig. 13b, it takes the real robot longer distances to make the sharp turning. This is also due to the larger ground resistance forces in reality.

3) Robustness to External disturbance:
We also test the FOC-PPOC-CPG controller's ability to keep tracking the desired target when the robot is disturbed by an external pushing force. Figure 14a and video "Disturbance recovery.mp4" shows an example trajectory of a disturbed goal-reaching task. It is observed from Fig. 14c that the FOC-PPOC-CPG controller reacts accordingly to its situation during the locomotion. When the deviation between the robot's head and the goal is relatively smaller before the disturbance (before 4.1s), the robot gently oscillates and adjusts its turning direction gradually towards the goal-direction. At around 9.04s, when the robot is pushed away from its desired direction, one can observe a clear redirection to the left-hand side of the robot's heading direction. The FOC-PPOC-CPG is able to adjust and make sharp turning to return to the correct direction and still reach the goal without wasting too much time on the recovery.

VII. CONCLUSION
This paper develops a bio-inspired controller for learning agile serpentine locomotion with a CPG net mimicking the central nervous system of natural snakes. The contribution of this paper is two-fold: First, we investigate the properties of the Matsuoka oscillator for achieving diverse locomotion skills in a soft snake robot. Second, we construct a PPOC-CPG net that uses a CPG net to actuate the soft snake robot, and a neural network to efficiently learn a closed-loop near-optimal control policy that utilizes different oscillation patterns in the CPG net. This learning-based control scheme shows promising results in goal-reaching and tracking tasks in soft snake robots.
This control scheme can be applicable to a range of biomimic motion control for soft robotic systems and may require the detailed design of the CPG network given insights from the corresponding biological systems. Our current research focuses on introducing sensory inputs into the CPG system, which enables reactive responses to contact forces with the external environment and generates an obstacle-aided locomotion controller for the soft snake robot. We will also investigate a scalable control design for extending the locomotion controller for a three-dimensional soft snake robot and a longer (more than 8 links) soft snake robot.

A. Describing function analysis of the Matsuoka Oscillator
According to Fourier theory, assume the vanilla state x(t) can be approximated by its main sinusoidal component, which is where r = d/A, r ∈ R is the ratio of bias to the amplitude of the signal. Then the output state z(t) = g(x(t)) = max (x(t), 0) can be written as: and The derivation of K(r) and L(r) are based on Fourier series analysis. Both K(r) and L(r) are constrained by −1 ≤ r ≤ 1 for x(t) to be non-negative in the period [−π, π]. The derivation of K(r) and L(r) can be found in Appendix A-B. Considering single period [− π ω , π ω ] for z(t), Given −1 ≤ r ≤ 1, when x(t) > 0, we have the remaining term (t) of the Fourier expansion simplified as . Based on the constraints, we can numerically calculate the bound of (t) for certain A and ω. For instance, if A = 1 and ω = 1, we have . Then by taking subtraction between flexor and extensor in eq. (1) and neglect small coupling terms from other primitive CPGs, we have which can be simplified to If x e i and x f i satisfy the perfect entrainment assumption, such that the amplitude A(x e i ) = A(x f i ) = A x , and the bias r(x e i ) = r(x f i ) = r x (both r x and A x are constants), and phase delay between x e i and x f i is π. In this case, we have z i = K(r x )x i . By setting K f = 1, eq. (A.6) can be further simplified to Next, an ordinary differential equation can be obtained by merging the two equations in eq. (A.7) as, When the system is harmonic, the coefficient of the firstorder derivative of eq. (A.8) becomes zero, then Specially, we use K n to represent K(r x ) for the harmonic Matsuoka oscillator.

B. Calculation of K(r) and L(r)
The general Fourier expansion of a continuous function g(x(t)) is: where a 0 = 1 π π −π g(x(t))dt g(x(t)) cos(t)dt g(x(t)) sin(t)dt.
Given the gate function g(x) = max(0, x), when the variable function x(t) = cos(t) + r, r ∈ R. In this case, the general function x(t) = cos(ωt) + r for any ω ∈ R has the same result as the special case with ω = 1. Such a case is trivial to prove, the reason for setting ω = 1 here is just for the space-saving purpose. Assume that the major component of g(x(t)) is also a cosine wave function, then it can be approximated as: g(x(t)) = 1 2 a 0 + a 1 cos(t).
When x(t) = cos(t) + r, both the bias a 0 and the amplitude a 1 become functions of r. Parameter r is the bias of the input wave signal. In this draft, we use K(r) to represent a 1 and L(r) to represent a 0 , which are calculated as follows: .
Specifically, according to the form of Fourier expansion, we have Where arccos r + arccos (−r) = π.
Specifically, in this case the coefficient of sine terms can be reduced to zero, This is because g(x(t)) is even and sin (ωt) is odd in this integral.

C. Amplitude Threshold of Transition from Free Oscillation to Forced Entrainment
According to Matsuoka's previous work on the frequency response of the Matsuoka oscillator [26], the free-response oscillation amplitude of the Matsuoka oscillator can be approximated as .
As long as τ r τ a ω 2 n − 1 > 0 is satisfied (in this paper the condition is satisfied, the readers may refer to the data in Appendix C to check the computation), function Q(U ) has its minimum value at U 0 = (τ r τ a ω 2 n − 1)(τ r τ a ω 2 − 1) + K 2 n τ 2 a a 2 ω 2 (τ r τ a ω 2 n − 1) + K 2 n τ 2 a a 2 ω 2 , (A. 21) where ω > ω n during the RL learning process. Then we have Since ψ i = z e i − z f i = g(x e i ) − g(x f i ) ≈ K n x i , we have ψ i = K n x i = K nxi + K n bias(x i ) = K nxi + K n 2u e − 1 (b − a)K n + 1 .

B. Proof of Proposition 2
Proof. If x e i and x f i do not satisfy the perfect entrainment assumption, then the amplitude A(x e i ) = A(x f i ), and bias r(x e i ) = r(x f i ). For simplification we denote A q := A(x q i ) and r q := r(x q i ) for q ∈ {e, f }. Instead of looking into the relation between x i and ψ i , we focus on the bias between the two states.

APPENDIX C DATA
This section includes the parameter configuration of the Matsuoka CPG network and the hyper parameter setting of domain randomization for the PPO algorithm.