Evolutionary Algorithm for the Theo Jansen Walking Mechanism

There is a Dutch artist/engineer who created a very elaborate walking mechanism. The working principle can be seen here:

http://www.strandbeest.com/beests_leg.php

The curious part is that he used a self-made Evolutionary Algorithm to calculate the ideal link lengths, which are described at the bottom of the page.

I created a Python script to visually analyze the ground-contact part of the cycle, which must have two requisites fulfilled:

1. Be as straight as possible, so as not to wobble up and down;
2. Have a speed as constant as possible, so as not to drag one foot against the other;

These two criteria would result in a "wheel like" effect, with the machine going linearly ahead without wasting kinetic energy.

The question is:

"Do you have any suggestion of a simple evolutionary iterative formula to optimize leg lengths (by inserting the correct mutations in the code below) so as to improve the walking path given the two criteria above?"

EDIT: some suggestions about the "fitting rule" for genome candidates:

• Take the "lower part" (ground contact) of the cycle, given that it corresponds to one third of crank revolution (mind the lower part might have a non-horizontal slope and still be linear);
• Apply linear regression on the point positions of this "ground contact" part;
• Calculate vertical variation from the linear regression (least squares?)
• Calculate speed variation by the difference between one point and the previous one, parallel to the regression line;
• (optional) plot graphs of these "error functions", possibly allowing to select mutants visually (boooring... ;o).

Here is my code, in Python + GTK, which gives some visual insight into the problem: (EDIT: now with parametrized magic numbers subject to mutation by mut's values)

# coding: utf-8

import pygtk
pygtk.require('2.0')
import gtk, cairo
from math import *

class Mechanism():
    def __init__(s):
        pass

    def assemble(s, angle):

        # magic numbers (unmutated)
        mu = [38, 7.8, 15, 50, 41.5, 39.3, 61.9, 55.8, 40.1, 39.4, 36.7, 65.7, 49]

        # mutations
        mut = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

        # mutated
        mn = [mu[n]+mut[n] for n in range(13)]

        s.A = Point(0,0)
        s.B = Point(-mn[0], -mn[1])

        s.C = fromPoint(s.A, mn[2], angle)
        s.ac = Line(s.A, s.C)

        s.D = linkage(s.C, mn[3], s.B, mn[4])
        s.cd = Line(s.C, s.D)
        s.bd = Line(s.B, s.D)

        s.E = linkage(s.B, mn[5], s.C, mn[6])
        s.be = Line(s.B, s.E)
        s.ce = Line(s.C, s.E)

        s.F = linkage(s.D, mn[7], s.B, mn[8])
        s.df = Line(s.D, s.F)
        s.bf = Line(s.B, s.F)

        s.G = linkage(s.F, mn[9], s.E, mn[10])
        s.fg = Line(s.F, s.G)
        s.eg = Line(s.E, s.G)

        s.H = linkage(s.G, mn[11], s.E, mn[12])
        s.gh = Line(s.G, s.H)
        s.EH = Line(s.E, s.H)


        return s.H


class Point:
    def __init__(self, x, y):
        self.x, self.y = float(x), float(y)
    def __str__(self):
        return "(%.2f, %.2f)" % (self.x, self.y)

class Line:
    def __init__(self, p1, p2):
        self.p1, self.p2 = p1, p2
    def length(self):
        return sqrt((p1.x-p2.x)**2 + (p1.y-p2.y)**2)

def fromPoint(point, distance, angle):
    angle = radians(angle)
    return Point(point.x + distance * cos(angle),
        point.y + distance * sin(angle))

def distance(p1, p2):
    return sqrt( (p1.x - p2.x)**2 + (p1.y - p2.y)**2 )

def ccw(p1, p2, px):
    """ Test if px is at the right side of the line p1p2 """
    ax, ay, bx, by = p1.x, p1.y, p2.x, p2.y
    cx, cy = px.x, px.y
    return (bx-ax)*(cy-ay)-(by-ay)*(cx-ax) < 0

def linkage(p1, l1, p2, l2):
    l1 = float(l1)
    l2 = float(l2)
    dx,dy = p2.x-p1.x, p2.y-p1.y
    d = sqrt(dx**2 + dy**2)                             # distance between the centers
开发者_JAVA百科    a = (l1**2 - l2**2 + d**2) / (2*d)                  # distance from first center to the radical line
    M = Point(p1.x + (dx * a/d), p1.y + (dy * a/d))     # intersection of centerline with radical line
    h = sqrt(l1**2 - a**2)                              # distance from the midline to any of the points
    rx,ry = -dy*(h/d), dx*(h/d)
    # There are two results, but only one (the correct side of the line) must be chosen
    R1 = Point(M.x + rx, M.y + ry)
    R2 = Point(M.x - rx, M.y - ry)
    test1 = ccw(p1, p2, R1)
    test2 = ccw(p1, p2, R2)
    if test1:
        return R1
    else:
        return R2




###############################33

mec = Mechanism()
stepcurve = [mec.assemble(p) for p in xrange(360)]

window=gtk.Window()
panel = gtk.VBox()
window.add(panel)
toppanel = gtk.HBox()
panel.pack_start(toppanel)

class Canvas(gtk.DrawingArea):
    def __init__(self):
        gtk.DrawingArea.__init__(self)
        self.connect("expose_event", self.expose)

    def expose(self, widget, event):
        cr = widget.window.cairo_create()
        rect = self.get_allocation()
        w = rect.width
        h = rect.height
        cr.translate(w*0.85, h*0.3)
        scale = 1
        cr.scale(scale, -scale)
        cr.set_line_width(1)

        def paintpoint(p):
            cr.arc(p.x, p.y, 1.2, 0, 2*pi)
            cr.set_source_rgb(1,1,1)
            cr.fill_preserve()
            cr.set_source_rgb(0,0,0)
            cr.stroke()

        def paintline(l):
            cr.move_to(l.p1.x, l.p1.y)
            cr.line_to(l.p2.x, l.p2.y)
            cr.stroke()

        for i in mec.__dict__:
            if mec.__dict__[i].__class__.__name__ == 'Line':
                paintline(mec.__dict__[i])

        for i in mec.__dict__:
            if mec.__dict__[i].__class__.__name__ == 'Point':
                paintpoint(mec.__dict__[i])

        cr.move_to(stepcurve[0].x, stepcurve[0].y)
        for p in stepcurve[1:]:
            cr.line_to(p.x, p.y)
        cr.close_path()
        cr.set_source_rgb(1,0,0)
        cr.set_line_join(cairo.LINE_JOIN_ROUND)
        cr.stroke()

class FootPath(gtk.DrawingArea):
    def __init__(self):
        gtk.DrawingArea.__init__(self)
        self.connect("expose_event", self.expose)

    def expose(self, widget, event):
        cr = widget.window.cairo_create()
        rect = self.get_allocation()
        w = rect.width
        h = rect.height

        cr.save()
        cr.translate(w/2, h/2)

        scale = 20
        cr.scale(scale, -scale)

        cr.translate(40,92)

        twocurves = stepcurve.extend(stepcurve)

        cstart = 305
        cr.set_source_rgb(0,0.5,0)
        for p in stepcurve[cstart:cstart+121]:
            cr.arc(p.x, p.y, 0.1, 0, 2*pi)
            cr.fill()

        cr.move_to(stepcurve[cstart].x, stepcurve[cstart].y)
        for p in stepcurve[cstart+1:cstart+121]:
            cr.line_to(p.x, p.y)
        cr.set_line_join(cairo.LINE_JOIN_ROUND)
        cr.restore()
        cr.set_source_rgb(1,0,0)
        cr.set_line_width(1)
        cr.stroke()




        cr.save()
        cr.translate(w/2, h/2)
        scale = 20
        cr.scale(scale, -scale)
        cr.translate(40,92)

        cr.move_to(stepcurve[cstart+120].x, stepcurve[cstart+120].y)
        for p in stepcurve[cstart+120+1:cstart+360+1]:
            cr.line_to(p.x, p.y)
        cr.restore()
        cr.set_source_rgb(0,0,1)
        cr.set_line_width(1)
        cr.stroke()



canvas = Canvas()
canvas.set_size_request(140,150)
toppanel.pack_start(canvas, False, False)

toppanel.pack_start(gtk.VSeparator(), False, False)

footpath = FootPath()
footpath.set_size_request(1000,-1)
toppanel.pack_start(footpath, True, True)


def changeangle(par):
    mec.assemble(par.get_value()-60)
    canvas.queue_draw()
angleadjust = gtk.Adjustment(value=0, lower=0, upper=360, step_incr=1)
angleScale = gtk.HScale(adjustment=angleadjust)
angleScale.set_value_pos(gtk.POS_LEFT)
angleScale.connect("value-changed", changeangle)
panel.pack_start(angleScale, False, False)


window.set_position(gtk.WIN_POS_CENTER)
window.show_all()
gtk.main()

Try the demo!

This is a fascinating question, though I think somewhat beyond the scope of Stack Overflow: it's not something that going to be solved in a few minutes, so I'll stick an outline here and update it if I make any progress. There are going to be three parts to any approach:

1. Scoring the footprint: does the linkage break? does the footprint have the right kind of shape? how flat is it? how smooth is the motion? does it spend enough time in the flat portion?

2. Searching for good values of the magic numbers. It's not clear that this has to be an evolutionary algorithm (though I can see why the idea of such an algorithm would appeal to Theo Jansen as it fits in with the animal metaphor in his art); perhaps other approaches like local search (hill climbing) or simulated annealing would be productive.

3. Searching for good configurations of the arms. This is where an evolutionary approach seems like it might be most worthwhile.

You can try out different magic numbers in my Javascript/canvas demo to see what kinds of motion you can get (CD = 55.4 is quite entertaining, for example). There's a whole mathematical theory of linkages, by the way, that connects the configuration spaces of linkages to topological manifolds.

I added some simple scoring to the demo. The ground score is the fraction of the cycle that the foot spends on the ground, which I take to be all points whose y-coordinate is within some tolerance of the lowest point. The drag score is the biggest difference between any two horizontal velocities while the foot is on the ground. (It's always negative, so that higher values = small differences in velocities = better.)

But here's where the difficulty comes in. In order to program any kind of search, I need to be able to combine these scores. But how do I balance them against each other? The Jansen magic numbers give me groundScore: 0.520; dragScore: -0.285. If I set AC=10, GH=65, EH=50, i get groundScore: 0.688; dragScore: -0.661. Nearly 70% of the time the foot is on the ground. But the take-off is draggy. Is it better or worse than Jansen's?

I think that getting actual engineering feedback in order to determine a good score is going to be the big problem here, not the actual search.