Locomotion is the kind of problem that animals have solved gracefully and effortlessly while robots lag, sometimes comically, behind. Everyone has seen bipedal robots fail miserably in the simplest terrains or when trying to use hand-feet coordination to, for example, open a door. It seems complex; computing the joint angles and positions with enough generality leads to a formidable control problem. The solution, natural evolution tells us, can be rather simple.

Simple agents, complex behaviours III: Oscillators [Robotics Tutorial]

Also Read: Simple agents, complex behaviours II: Coincidences [Robotics Tutorial]

The control architecture of relatively high-order animals possesses a hierarchical organization that relies on two main aspects at its lowest level. First, as we have emphasized in this series, efficient body design makes it possible to generate locomotion without any actuator, motor or muscle just by exploiting the physics of rigid bodies. In this case, the energy source is gravity; we will not work out the details in this article, but the reader is encouraged to watch some videos to whet their intellectual appetite.

Second, when an energy source is required to do something other than gracefully descending a hill, animals count on Central Pattern Generators (CPG); they use oscillations. Consider the walk or the Toelt of a horse, its trot or its pace; all are different oscillation modes. They are rhythmical patterns that require minimal upstream control in a constant terrain. Even more, the fact that we can categorically classify them indicates that there is not much in the middle; they are all abrupt changes in the dynamics of the underlying system; they are bifurcations.

 

The World of Non-linear Oscillations

The simplest model for oscillatory phenomena is the sine function; in its crests and troughs, it leads us on an endless trip around the circle. Even though we can use it as a building block in our locomotion, it falls short in the number of different behaviours it can display. We need a more powerful kind of oscillation, a non-linear oscillation. While the crests and troughs in the sinusoidal world are symmetrical, non-linear oscillators are asymmetrical and display different time scales, among other exciting features.

Figure 1 shows examples of different modes of a non-linear oscillator. The transition between them is a bifurcation: a qualitative change in behaviour due to changes in the model's parameters. Let's explore these concepts by building our own CPG-based locomotion system.

 

Different Energy States

When studying non-linear oscillators, it is useful to start at an energy landscape. In the real and mathematical worlds, many systems enjoy low-energy states and chase them as long as it is possible. This gives us an intuitive picture of what a non-linear oscillation is. Consider below figure:-

We locate the system's state in the abscissas and its energy in the ordinates; the blue curve is the energy landscape. If we initialize the system at point A, it will quickly roll until a minimum is reached at B. However, suppose we deform the curve so that B disappears. In that case, the system will transition to state C. These are two primary components of our oscillation: fast rolling dynamics that settles on local minima and a slow deformation that extinguishes those minima.

We made it this far without equations! Quite a feat. But we can only get this far without them; we need to provide an evolution law for our fast and slow variables.

 

Slow and fast motions: The Van der Pol oscillator

First, we will use the following double well potential as our energy landscape:

which we can plot using the following code:-

Our dynamical system, as we already hinted, will have two variables. The first one will be in charge of the "rolling down" we mentioned, mathematically this means that it will follow the gradient downwards:-

This is the fast variable, as the rolling happens rather fast relative to the deformation. We invite the reader to plot different values of the parameter "y". Notice that the shape changes as we vary the parameter! Parameter "y" is then a good candidate for our slow variable, which we evolve using a simple linear decay proportional to “x”:

And that's it; we have our non-linear oscillator with fast-slow dynamics. This particular model is called the Van der Pol oscillator. Let's write it down as a system of Ordinary Differential Equations:-

The parameter ε < 1 is very small, which guarantees the system's different time scales. Let's see how it works.

 

The Dynamics

The first step in the analysis is to find the nullclines of both variables. The nullclines are the curves at which evolution stops for each of our variables, that is, when the right-hand side is equal to zero. The first nullcline is called the critical manifold and, if you think about it, shows precisely the location of the minima of the potential as a function of the parameter "a".

Here is the code to display the below curve:-

import numpy as np
import matplotlib.pyplot as plt

b = 0.0
y1 = lambda x: -x + x**3
x2 = lambda y: b*np.ones_like(x)

x = np.linspace(-1, 1, 100)
y = np.linspace(-1, 1, 100)
plt.plot(x, y1(x), label = 'x-nullcline')
plt.plot(x2(y), y, label = 'y-nullcline')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

This curve has three sections shown in the figure. The important thing to remember is that the system really doesn't like to be in the middle branch; let's use Euler's method to find a solution and see what we mean:-

The result is fascinating:-

When the trajectory is away from the critical manifold, it quickly jumps towards it, hugs it and stays there as long as possible. But all good things come to an end, and when it reaches the fringes of the evil middle branch, it jumps, hoping to reach a more comfortable place to continue its love affair. This completes the oscillation.

Plotting the variables in the time domain reveals the final form of the wave. Notice the plateaus and fast transitions. Can you anticipate why they are important in a walking agent? An interesting aspect of the dynamics is the phenomenon of hysteresis; the jumping point of the trajectory going up is not the same as the one for the trajectory going down!

Here is the complete code:-

import numpy as np
import matplotlib.pyplot as plt

T = 1000
b = 0.0
h = 0.01
epsilon = 0.05
# Evolution laws
f = lambda x, y: (y + x - x**3)/epsilon
g = lambda x, y: b - x
# Nullclines
y1 = lambda x: -x + x**3
x2 = lambda y: b*np.ones_like(y)

x0 = 0.0
y0 = 1.0
x = np.zeros(T)
y = np.zeros(T)
x[0] = x0
y[0] = y0
t = np.zeros(T)

for i in range(T-1):
    x[i + 1] = x[i] + h*f(x[i], y[i])
    y[i + 1] = y[i] + h*g(x[i], y[i])
    t[i + 1] = t[i] + h

xx = np.linspace(-1, 1, 100)
yy = np.linspace(-1, 1, 100)
plt.plot(xx, y1(xx), label = 'x-nullcline')
plt.plot(x2(yy), yy, label = 'y-nullcline')
plt.plot(x, y, 'k', linewidth = 2.0)
plt.ylabel('y')
plt.xlabel('x')
plt.show()

 

Different Gaits

So far, we only talked about the first nullcline; remember that those curves are the points at which either of the two variables stops evolving. When we add the second nullcline, we can graphically identify the system's fixed points. With it, we can gain some understanding of how the system transitions between different behaviours; those transitions are what we call bifurcations.

Notice that in the left part of the figure, the two nullclines intercept in the middle branch of the first nullcline and that branch, we established, is unstable, so this point is unstable as well, meaning every trajectory wants to stay away from it (we can see this by noting that looking at the vector field nearby that point but that is beyond the goal of this article). In the sequence of figures to the right, we can see the effect of changing the parameter "b"; the point of intersection moves to the right until the oscillation disappears when such a point is located in one of the stable branches. This is called a Hopf bifurcation.

Here is the complete code to display above figure:-

import numpy as np
import matplotlib.pyplot as plt

T = 1000
h = 0.01
epsilon = 0.05

def integrate( ax, b ):
    # Evolution laws
    f = lambda x, y: (y + x - x**3)/epsilon
    g = lambda x, y: b - x
    # Nullclines
    y1 = lambda x: -x + x**3
    x2 = lambda y: b*np.ones_like(y)

    x0 = 0.0
    y0 = 1.0
    x = np.zeros(T)
    y = np.zeros(T)
    x[0] = x0
    y[0] = y0
    t = np.zeros(T)

    for i in range(T-1):
        x[i + 1] = x[i] + h*f(x[i], y[i])
        y[i + 1] = y[i] + h*g(x[i], y[i])
        t[i + 1] = t[i] + h

    xx = np.linspace(-1, 1, 100)
    yy = np.linspace(-1, 1, 100)
    ax.plot(xx, y1(xx), label = 'x-nullcline')
    ax.plot(x2(yy), yy, label = 'y-nullcline')
    ax.plot(x, y, 'k', linewidth = 2.0)
    ax.set_ylabel('y')
    ax.set_xlabel('x')

fig, ax = plt.subplots(1,3)

integrate(ax[0], 0.0)
integrate(ax[1], 0.579)
integrate(ax[2], 0.9)
plt.show()

 

A Walking Agent

As we said before, CPGs can be used to implement locomotion in artificial agents similarly to natural agents. With our non-linear oscillator, we can readily implement two-legged locomotion by connecting its output to some simple muscles. Here is our setup. First, let's create a class to encapsulate the CPG:-

Then, we create some simple links that represent the agent's extremities:-

Finally, we iterate the model initializing the CPGs at different phases with respect to each other:-

As we can see, we can implement at least two different gaits in our agent.

The first one (top) is a regular pace walk in a constant terrain. This is achieved when the oscillator has a linear behaviour (a = 0.1, b = 0.4), as in a sinusoidal wave, providing the environment that allows for a constant grip on the actuators. The second one (bottom) is a slower walk that maximizes the contact time with the ground (a = 1.0, b = 0.4); this is useful in more slippery terrains.

Here is the complete code:-

# Implementation of a simple walking pattern
import numpy as np
import matplotlib.pyplot as plt

class CPG:

   def __init__(self, x0, y0, b):
       self.b = b
       a = 1.0
       self.h = 0.01
       self.epsilon = 0.05
       self.f = lambda x, y: (a*x - x**3 + y)/self.epsilon
       self.g = lambda x, y: self.b - x
       self.t = 0.0
       # State
       self.x = x0
       self.y = y0

   def step( self ):
       dx = self.f(self.x, self.y)
       dy = self.g(self.x, self.y)
       self.x = self.x + self.h*dx
       self.y = self.y + self.h*dy
       self.t = self.t + self.h

       return self.x


class Link:

   def __init__(self, x0, y0, cpg, theta_min, theta_max, l = 0.5):
       self.x0 = x0
       self.y0 = y0
       self.theta = 0.0
       self.theta0 = -np.pi/2.0
       self.theta_min = theta_min
       self.theta_max = theta_max
       self.l = l
       self.cpg = cpg

   def step( self ):
       r = self.cpg.step()
       xu = 1.0
       xl = -1.0
       m = (self.theta_max - self.theta_min)/(xu - xl)
       self.theta = m*(r - xl) + self.theta_min


   def draw(self, ax):
       x = self.l*np.cos(self.theta0 + self.theta)
       y = self.l*np.sin(self.theta0 + self.theta)
       ax.plot([self.x0, self.x0 + x], [self.y0, self.y0 + y], 'k')

       links = [Link(1.0, 1.0, CPG(1.0, 0.0, 0.4), -0.2, 0.0), 
       Link(1.0, 1.0, CPG(-1.0, 0.0, 0.4), -0.2, 0.0), 
       Link(1.0, 1.5, CPG(1.0, 0.0, 0.2), -0.1, 0.1, 0.4), 
       Link(1.0, 1.5, CPG(-1.0, 0.0, 0.2), -0.1, 0.1, 0.4)]

       fig, ax = plt.subplots(1,1)
       plt.show(block = False)
       x = -0.8

       for i in range(1000):
       # ax.cla()

       x += 0.05

       for link1 in links:
       link1.x0 = x 
       link1.step()
       if np.mod(i, 10) == 0:
           link1.draw(ax)

       if np.mod(i, 10) == 0:
           ax.plot([x, x], [1.0, 1.5], 'k')
           plt.pause(0.01)

       plt.axis([0, 20, 0, 2])

plt.show()

 

Conclusions

In this article we have implemented a simple CPG for locomotion. This is a toy example that is overly simplified to show the fundamentals of non-linear oscillators. In a real control architecture, the parameters of such a dynamical system are controlled by high level mechanisms that modulate the bifurcations with a sensory-motor feedback loop.

 

References

1. https://en.wikipedia.org/wiki/Passive_dynamics
2. Katz, P. S. (2016). Evolution of central pattern generators and rhythmic behaviours. Philosophical Transactions of the Royal Society B: Biological Sciences, 371(1685), 20150057.
3. Nassour, J., Hénaff, P., Benouezdou, F., & Cheng, G. (2014). Multi-layered multi-pattern CPG for adaptive locomotion of humanoid robots. Biological cybernetics, 108(3), 291-303.