Minimum time value iteration solution of the double integrator

I've spent the last week of spare time working of a value iteration problem for the underactuated robotics subject I've been working on. In particular this has involved finding the optimal policy of how to get the undamped double integrator to the origin in minimum time (github gist of code). A MATLAB / Octave solution is available on the OCW site, so I don't feel too bad publishing this here.

This has been something a bit different for me, considering how little vectorised code I've written in the past, but has already had a big impact on some things I'm doing for a number of projects. The performance benefit is substantial, but it makes code hard to read, write, debug and test.

Python simulation of the simple pendulum

Following on from the previous post, I've created some simulations of the simple pendulum without damping or control input. This work was designed more to be a test of my understanding of how to use the scipy integrator module, and to make sure the matplotlib plotting code worked correctly.

I'm pretty pleased with the results, and you can see a plot of the homoclinic orbits of pendulum starting in different initial states, and some videos of what that means, if you don't naturally think in phase space.

Python simulation of the Van der Pol Oscillator

One of the things I'm doing at the moment is watching an excellent series of MIT OCW lectures by Assoc. Prof. Russ Tedrake on underactuated robotics. I've learned a lot watching the videos, which pick up at a point where some parts of my doctoral studies finished.

Just to pick something small that is at least tangentially related to the related to the course, I've attached a python snippet of code for generating the Van der Pol Oscillator. This system is significant in that it exhibits limit cycles, which are an important tool for reasoning about walking robots

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

mu = 0.2

def van_der_pol_oscillator_deriv(x, t):
    nx0 = x[1]
    nx1 = -mu * (x[0] ** 2.0 - 1.0) * x[1] - x[0]
    res = np.array([nx0, nx1])
    return res

ts = np.linspace(0.0, 50.0, 500)

xs = odeint(van_der_pol_oscillator_deriv, [0.2, 0.2], ts)
plt.plot(xs[:,0], xs[:,1])
xs = odeint(van_der_pol_oscillator_deriv, [-3.0, -3.0], ts)
plt.plot(xs[:,0], xs[:,1])
xs = odeint(van_der_pol_oscillator_deriv, [4.0, 4.0], ts)
plt.plot(xs[:,0], xs[:,1])
plt.gca().set_aspect('equal')
plt.savefig('vanderpol_oscillator.png')
plt.show() 

Which when simulated looks something like this: