Hanayama O'Gear Puzzle

Declan showed me a brain teaser his friend sent him. Apparently, it's known as the Hanayama O'Gear Puzzle. It consists of a bronze colored metal cube with a captive gear which can roll from face to face on the cube. The gear can also rotate between two directions on each face. Some edges of the cube are made so that the gear can't roll over them. One tooth of the gear and one face of the cube have a special cut-out. Only when the special tooth is engaged with the special face can the gear be freed. The object of the puzzle is to manipulate the gear into this orientation. Here's a video of the puzzle in action.

I tried solving the puzzle for about fifteen minutes before deciding to use a computer. My basic strategy was to treat the problem as a graph search, where the nodes of the graph are the different states the cube and gear can be in. At any time, you can only do three things: rotate, roll forward, and roll backward. I figured that the limited set of moves and possible positions would prevent the search space from getting too big.

I marked up the cube with stickers, naming each face and defining its 'north' direction. I also marked the gear by labeling each tooth and setting a forward rotation direction.

I don't know Haskell very well, so I thought using it would be a fun challenge. Plus, I had Haskell guru Phil to ask for help and advice. A lot of the code is just defining the connectivity of the puzzle. I used a list to keep track of the links between different faces and a pattern-matched function to keep track of the possible rotations. I defined north on each face so that these rotations were always the same. A north-facing gear could only rotate to the east, and a south-facing one could only rotate west.

Beyond the trickiness of defining directions and connectivity, I just implemented a straightforward breadth-first search. The program creates a list of all states accessible from the starting state. Because the list is generated with successive expansions of the search frontier, it is naturally sorted by search depth. I take advantage of Haskell's lazy evaluation, creating an infinite list of accessible states but only ever calculating a few of them. Code below.

import Data.Maybe (catMaybes)
import Data.List (find)
import Data.Tuple (swap)

data Face = Top | Bot | Lef | Ri | Front | Back deriving (Show, Eq)
data Dir = North | South | East | West deriving (Show, Eq)

data GearState = GearState{ face :: Face, dir :: Dir, tooth :: Int, history :: [(Face, Dir)] } deriving Show

-- a list defining connections between faces
forward_connections = 
    [ ((Top, East),(Ri, South))
    , ((Top, West),(Lef, South))

    , ((Bot, North),(Front, West))
    , ((Bot, South),(Back, East))
    , ((Bot, West),(Ri, North)) 

    , ((Front, South),(Bot, North)) 
    , ((Front, West),(Back, West)) 

    , ((Back, South),(Bot, North)) 
    , ((Back, East),(Ri, East)) 
    , ((Back, West),(Lef, West)) 

    , ((Lef, North),(Top, East)) 
    , ((Lef, East),(Front, East))
    , ((Lef, West),(Back, West)) 
    , ((Ri, North),(Top, West)) 
    , ((Ri, South),(Bot, East)) 
    , ((Ri, East),(Back, East))  

-- mapping the reverse tuple function over the connections gives the connections for reverse gear moves
backward_connections = map swap forward_connections

-- a function which returns a gear moved forward or backward
-- takes a list of connections between faces and a function for finding the next tooth
-- forward_connections and incr for forward moves
-- backward_connections and decr for backward moves
move :: GearState -> [((Face, Dir),(Face, Dir))] -> (Int -> Int) -> Maybe GearState
move (GearState face dir tooth history) connections nextTooth = 
    case (lookup (face, dir) connections) of -- looks for a connection
        Nothing -> Nothing -- if no forward connection exists from the current state
        (Just (new_face, new_dir)) -> Just (GearState new_face new_dir (nextTooth tooth) ((new_face,new_dir):history)) -- returns a new state, incrementing the tooth and adding to the history

-- a function which applies a turn to the gear
turn :: GearState -> GearState
turn (GearState face North tooth history) = GearState face East tooth ((face,North):history)
turn (GearState face East tooth history) = GearState face North tooth ((face,East):history)

turn (GearState face South tooth history) = GearState face West tooth ((face,South):history)
turn (GearState face West tooth history) = GearState face South tooth ((face,West):history)

-- some mod 5 math for the teeth
mod5 :: Int -> Int
mod5 = (flip mod) 5

incr :: Int -> Int
incr x = mod5 (x+1)

decr :: Int -> Int
decr x = mod5 (x-1)

-- determines if the gear can leave the cube
isWinningState :: GearState -> Bool
isWinningState (GearState Top North 0 _) = True
isWinningState _ = False

-- creates a list of all the states accessible from the current state
nextStates :: GearState -> [GearState]
nextStates state = catMaybes [Just (turn state), (move state forward_connections incr), (move state backward_connections decr)] -- removes Nothings and strips the values out of Just

-- creates a new search frontier by replacing each state in the search frontier with all the states accessible from that state and flattening
newFrontier :: [GearState] -> [GearState]
newFrontier frontier = (concatMap nextStates frontier)

-- creates a list of all future frontiers
reachableStates :: GearState -> [GearState]
reachableStates start = concat (iterate newFrontier [start])

-- finds the winning state
solve :: GearState -> Maybe [(Face, Dir)]
solve start = fmap history (find isWinningState (reachableStates start))

You can run it with a command like, solve (GearState Top South 4 []). This code, and a different Python solver I wrote are available in my Github.


Couch Table

Matt and I needed a quick project for Christmas time. We came up with a small table to fit over the arm of a couch.

It's an original design with a u-shaped base that allows the table to slide forward and out of the way of the arm rest by wrapping around the leg of the couch. We made it out of red oak with two small walnut accents. Plans below.

First we built the base. We ripped 2" nominal thick red oak into 1.5" strips and chopped them to length. We cut shallow angles into the parallel parts of the base using a band saw and hand plane. We used hand tools to cut mitered bridle joint corners. These were the toughest part of the project, and we spent a lot of time fussing over the precise fit.

We made the top of the table by chopping a piece of 1x12" red oak. We beveled three edges using a hand plane and cut the through mortices.

Next, we made the vertical legs of the table. We ripped 1.25" strips of 1" nominal red oak and chopped them to length. We cut through tenons for the top and bottom with 1/8" shoulders.

We decided to try wedged tenons for a nice look on the top of the table. By eye, we cut some slim wedges out of walnut on the band saw. We also cut slots into the top tenon and flaired out the corresponding mortices slightly to accomodate the extra width from the wedges.

We glued everything together and gently hammered the wedges in.

Here's what they looked like after we sawed off the excess.

Finally, we planed all of the joint surfaces flat and sanded down to 200 grit. For a finish, we decided to use danish oil. I highly recommend it. At least for red oak, it make the grain stand out and gives the whole thing a wonderful color. We also applied a layer of finishing wax.

This table was fun to build and came out looking good. It doesn't feel too unstable, either. With some simple length adjustments, it would be easy to fit this design to different couches and chairs.


Wedge Table

This year over Thanksgiving, Matt and I made an interesting table out of red oak and walnut. It was inspired by a clipping from a catalog. This was our hardest project so far, and I think it came out well.

The original design is by Eben Blaney. Thanks Eben! We worked from the picture above and a few others, and tried to roughly match his table.

We started by modeling with Onshape.

The design ended up being pretty complicated. We picked a 45° angle between the long leg and the top. Then, to make things easy and reduce the chance of cutting angles wrong, we made the central triangle of the table isosceles. That meant that every angle on the table was either 45° or 67.5°. I think this was a good idea.

First, we chopped all the wood to length. Two 38" lengths of 1x4" for the double legs, a 53.5" length of 1x8" for the single leg, and a 48" length of 1x12" for the top. We also cut some smaller pieces of walnut for the stretchers and the peg.

Next, we bent a thin piece of wood to mark the curves on the double legs and carefully cut the excess away with a jig saw. The jig saw blade wanders around and doesn't make edges that are perfectly square. We fixed this a bit using sand paper. I think a band saw would have been better for this, but we didn't have one.

Using a pull saw and chisels, we cut diagonal dados into the double legs where the single leg would cross it.

Cutting the angled hole in the top was tricky. We made a jig and roughed out this hole with a hand drill, then finished the job with chisels.

We cut the shoulders off of the single leg, fitted it into the hole in the top, marked the hole for the pin, and cut it using the same technique.

We finished the double legs by cutting the tenons for the stretchers. Normal tenons for the top stretcher and through tenons for the bottom one. We made the stretchers themselves out of walnut.

Next, we glued up the legs. Because the top was not glued to the legs, we added a pin between the top and the upper stretcher to keep the top from sliding around. We cut an angled peg for top by eye out of walnut.

Finally, we sanded everything and applied finishing wax.

Came out pretty good, I think.


Solving the Brachistochrone with PyTorch

The brachistochrone problem is simple to explain. Imagine a bead which is free to slide on a wire. What shape of wire will get the bead from point A to point B in the shortest time? If you're wondering about the name, "brachistochrone" apparently means "shortest time" in Ancient Greek. Makes sense. The image below shows the optimal path in red, along with some slower paths in blue.

This was inspired by Declan, who was using this problem to learn about genetic algorithms. To quote him, "The genetic optimization algorithm may be the greatest concept to ever emerge from the mind of man." A bit dramatic, I think.

The problem can be solved with calculus, but I decided to approach this problem with gradient descent using PyTorch. It seemed like a natural solution, and a good way to learn about Torch. The image below shows my optimization approaching the shape of the known cycloid solution in blue.

I divided the wire into segments with heights \(y_i\) and calculated the speed at each segment using the height difference from the starting point. Then I found the length of the segments using the height differences between points. I calculated the time to reach each point, \(t_i\), by integrating the segment length divided by the speeds. Finally, I used automatic differentiation to find the derivatives \(\frac{dt_i}{dy_j}\) and used gradient descent to minimize the time to reach the final point, \(t_{last}\).

import torch
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize

starting_y = 0.25 # y position of the beginning of the wire, end always at 0
xs = torch.linspace(0,1,25) # evenly spaced points in x, from 0 to 1
ys = torch.linspace(starting_y,0,25, requires_grad=True) # start with a straight line connecting start to end

# the physics
def times(ys): # calculate the time to get to each point
	vs = torch.sqrt(starting_y - ys) # find velocity based on height
	dys = ys[1:] - ys[:-1] # find y difference
	dx = xs[1]-xs[0] # find x  difference
	lengths = torch.sqrt(dx.pow(2) + dys.pow(2)) # calculate arc length
	return torch.cumsum((2./(vs[:-1]+vs[1:])) * lengths, dim=0) # integrate to find time to each point using midpoint velocity

fig = plt.figure()
ax = fig.add_subplot(111)

# solving for the cycloid solution
def func(theta):
	return (1.-np.cos(theta)) / (theta - np.sin(theta)) - .25

theta = optimize.brentq(func, 0.01, 2*np.pi)

# parameterized cycloid solution
r = 1 / (theta - np.sin(theta))
ts = np.linspace(0,2 * np.pi,100)
sol_xs = r * (ts-np.sin(ts))
sol_ys = starting_y - r * (1-np.cos(ts))

ax.plot(sol_xs, sol_ys)

li, = ax.plot(xs.numpy(), ys.detach().numpy())


# a pause

# gradient descent
lr = 0.005 # learning rate
for i in range(2000):
	ts = times(ys) # get the times to get to each point
	t = ts[-1] # the time to get to the final point
	t.backward() # back propagate to find gradients
	grads = ys.grad # hold on to the gradient

	# update plots

	# update wire heights
	with torch.no_grad():
		ys[1:-1] -= lr * grads[1:-1]
	ys.grad.data.zero_() # zero gradients

ZIPY — A Homebrew Inverted Pendulum and Control System

Phil and I have been working on this one for a long time. We're calling it ZIPY: Zippy Inverted Pendulum Yakshave. We finally succeeded in controlling an inverted pendulum. It swings up and balances. Would ya look at that! Also check out Phil's posts on the project.

The inverted pendulum or cart pole is a classic problem in control theory. It's in OpenAI Gym of course, but we wanted to see it work in real life, not some lame simulation. We were inspired by Russ Tedrake's class on underactuated robotics and by the popularity of the problem for reinforcement learning.

A Real-Life Cart Pole

It took a few iterations, but we eventually found a system that works well. Our cart is 3D printed PLA driven by a DC motor via a toothed belt. The pole itself is a paint stirrer. One of the longer type, about 24". A rotary encoder opposite the motor acts as a pulley for the belt and allows us to track the motion of the cart, while a second rotary encoder on the cart is a pivot for the pole and measures its angle. The motor is controlled by a 32 amp Sabertooth motor controller. It's overkill, and pretty expensive at about $120, but we already had it for another project. We monitored the encoders with an Arduino. The foundation of the system is a piece of extruded aluminum rail called V-Slot, on which the cart slides and the motor and encoder are mounted. Our rail is 1.5 m long, from a company called SMW3D.

We used a 3D printer to make a PLA cart with mounting points for a rotary encoder and grooves to lock in a toothed belt. We designed all of our parts in Onshape. It's my favorite.

We found that if you slide the cart onto the V-Slot and dunk the two in boiling water for a moment, you a good fit for sliding with no fancy bearings. A touch of 3-in-1 oil seems to help as well. We also printed some pieces to mount the motor and rotary encoder to the rail. All of the designs can be found in our Onshape project.

The Software

All the software can be found in our Git repositiory. The rotary encoders were monitored by an Arduino Uno. Rotary encoders produce pulses on their two outputs as the shaft rotates. The pulses occur at regular intervals as the shaft turns, and the relative timing of the pulses give the direction of rotation. Driving them is just a matter of counting and interpreting these pulses. Normal operation requires two interrupt pins per encoder, but you can get away with one if you're willing to accept half the resolution. That's what we did because the Uno only has two interrupt pins. We still got 600 positions per revolution of resolution. Good enough. The microcontroller counts these pulses and keeps track of position by dead reconing. The program sends the current encoder positions whenever it receives a query over serial.

On the other side, a Python program contains all the smarts. It commands the motor controller over serial and can poll the Arduino whenever it needs to know the state of the system. We used two different controllers, a controller to pump energy into the system and swing the cart into the upright position, and an linear quadratic regulator to balance it upright.

The LQR Controller

The balancing controller is called an linear quadratic regulator or LQR controller. Phil wrote a more detailed explanation of the math involved, but here's an overview. It's an optimal controller for a linear system in a state, $$x = \begin{bmatrix}q\\\dot{q}\\\theta\\\dot{\theta}\end{bmatrix},$$ with control input, $$u = \begin{bmatrix}f\end{bmatrix}$$ equations of motion, $$\dot{x} = Ax + Bu,$$ and a quadratic cost function, $$c = x^T Q x + u^T R u,$$ Where \(q\) is the position of the cart, \(\theta\) is the angle of the pole counter-clockwise from vertical, \(f\) is the force on the cart, \(A\) and \(B\) are matrices that define the dynamics, and \(Q\) and \(R\) are matrices that define the cost.

To obtain the LQR controller, we first found the equations of motion. We assumed that our pole was light enough that it didn't significantly affect the cart's motion, reducing the free parameters to just \(\theta\) and \(\dot{\theta}\). The only forces on the pole are due to gravity and the motor. The torque on the pole from gravity is \(\tau = \frac{mgl}{2} \sin{\theta}\), so the angular acceleration is \(\ddot{\theta} = \frac{\tau}{I} = \frac{m g l}{2I}\sin{\theta} \approx \frac{m g l}{2I}\theta\). The torque on the pole from the motor is identical to that of gravity, but turned 90 degrees in \(\theta\), so \(\ddot{\theta} = \frac{m \ddot{x} l}{2I}\cos{\theta}\ \approx \frac{m \ddot{x} l}{2I}\). The effects of gravity end up in \(A\), $$A = \begin{bmatrix}0 & 1 & 0 & 0\\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & \frac{m g l}{2I} & 0\end{bmatrix},$$ and those of the cart acceleration end up in \(B\), $$B = \begin{bmatrix} 0 & 1 & 0 & \frac{m l}{2I}\end{bmatrix}.$$ The controller \(K\) allows you to find the optimal control input for a state \(x\), $$ u = -K x.$$ You can find \(K\) by solving a nasty equation for \(P\), $$A^TP + PA - PBR^{-1}B^TP + Q = 0,$$ and plugging into, $$K= R^{-1} B^T P.$$ Luckily, there is a SciPy function that does it for you.

The Swing-Up Controller

A clever trick allows us to pump energy into the system when the energy is below the target energy (the energy of a motionless pendulum, standing straight up), and suck energy out when the energy is above the target. $$f = (E(x)-E_{target}) \dot{\theta} \cos(\theta)$$ Imagine the cart is below the energy threshold and moving counter-clockwise. You want to apply a counter-clockwise torque to increase the rotation speed. If the pole is above the horizontal, that means pushing the cart right with a positive force. If it is below the horizontal, the cart should be pushed left. This flip is captured by the \(\cos{\theta}\). Flip the sign of the diffence from target energy or the current direction of rotation, and you also want to flip the input. So overall, this equation always provides the right direction for the force. Pretty neat!


What didn't work

We tried monitoring the position of the cart and pole using an OpenCV vision system. We put large colorful stickers at the top and bottom of the pole, and looked for areas in the images with those colors. You can see in our images that we had all kinds of junk in frame along with the pendulum. By choosing the right colors and removing prolematic objects, we were able to make this work. It was very annoying, though. As light changed in the room, we had to recalibrate the system often.

There was a bigger problem, though. We discovered that OpenCV (or possibly the drivers of the webcams we tried) introduced a five frame buffer, which caused a five frame delay. We never figured out a way around this and we never succeeded in controlling the system with this delay.