Ramp.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:48325435f290ef432fc24801de0335ba9e6eb4635fbd2532838155deb6d1e652"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## 3D printable candidate solutions to the brachistochrone problem\n",
      "\n",
      "Consider the problem of finding a curve that connects two points A and B,\n",
      "such that an object rolling down the curve from A to B takes the least amount of time to\n",
      "make the trip\n",
      "\n",
      "![Alt text](static/Brachistochrone.png)\n",
      "\n",
      "This is called the [Brachistochrone problem](http://en.wikipedia.org/wiki/Brachistochrone_curve), a classic problem in the [Calculus of Variations](http://en.wikipedia.org/wiki/Calculus_of_variations).\n",
      "\n",
      "This notebook contains functions which will let you set up your own version of the problem, sketch out your own solution using your mouse, then generate a 3D-printable ramp to roll marbles down in order to test your design.\n",
      "\n",
      "Source code for all func\n",
      "\n",
      "### Requirements\n",
      "\n",
      "The functions of this notebook require:\n",
      "\n",
      "- [repository of source files associated with this notebook](https://github.com/thearn/printable_ramps)\n",
      "- ipython\n",
      "- numpy\n",
      "- scipy\n",
      "- matplotlib\n",
      "- solidpython\n",
      "- OpenCV\n",
      "- OpenSCAD\n",
      "\n",
      "The python libraries can be installed by running `pip install -r requirements.txt` within the directory that contains this notebook file. \n",
      "\n",
      "The computer vision library OpenCV can be installed [by following their installation instructions](http://docs.opencv.org/doc/tutorials/introduction/table_of_content_introduction/table_of_content_introduction.html). OpenCV can be tricky to install. I'd recommend using a system package manager like apt-get or homebrew if possible. On Windows, the easiest way to get started with OpenCV and Python is to use either [Python(x,y)](https://code.google.com/p/pythonxy/) or the [Anaconda Python Distribution](https://store.continuum.io/cshop/anaconda/)\n",
      "\n",
      "The CAD program OpenSCAD can be [downloaded from here](http://www.openscad.org/downloads.html). This notebook involves calling OpenSCAD from the command line. To make this work, the path to the OpenSCAD binary will need to be given in the next cell.\n",
      "\n",
      "## Getting started\n",
      "\n",
      "First, let's make all necessary imports and set up the conditions (the location of the start and end points A and B) for the problem. Feel free to modify these points if you'd like.\n",
      "\n",
      "The path to the OpenSCAD (including the name of the binary) needs to be given below. The default location given is the path on my Mac, you'll need to change this if necessary. If OpenSCAD is already in your system's path, the binary name should be all that is needed here."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from lib.ramp import ramp, render_stl\n",
      "from lib.interface import get_sketch, smooth_plot, show_stl, Animate_frame\n",
      "from lib.fit import fit_multiple, rescale\n",
      "from lib.cycloid import cycloid, total_time\n",
      "\n",
      "from IPython.html.widgets import interactive, FloatSliderWidget\n",
      "from IPython.display import display\n",
      "from IPython.display import HTML, Image\n",
      "import os\n",
      "\n",
      "openscad_path = \"/Applications/OpenSCAD.app/Contents/MacOS/OpenSCAD\"\n",
      "\n",
      "# tests the given path\n",
      "if not os.path.exists(openscad_path):\n",
      "    raise Exception(\"OpenSCAD path is incorrect!\")\n",
      "\n",
      "A = [0, 95.0] # Start point\n",
      "B = [210, 55.0] # End point\n",
      "\n",
      "n_points = 250 # resolution"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Run the next cell, and draw a curve with your mouse between the two dots shown in the window.\n",
      "\n",
      "(Optional) Use the widget below to smooth out your drawn curve after it has been collected. If sigma = 0 (default), no smoothing will be performed."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x,y = get_sketch(A, B, n_points)\n",
      "\n",
      "y_smooth = y\n",
      "def f(sigma):\n",
      "    global y_smooth\n",
      "    y_smooth = smooth_plot(x, -y, sigma, A, B)\n",
      "\n",
      "sigma_slider = FloatSliderWidget(min=0, max=15, step=0.3, value=0)\n",
      "w=interactive(f,sigma=sigma_slider)\n",
      "display(w)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now collect some more curves that connect A to B to compare against, plot them, and compute a prediction of descent time for each curve.\n",
      "\n",
      "Move the slider to simulate a small object rolling down each curve over time."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Generate an inverted cycloid between A and B with a given number of sample points\n",
      "x, y0 = cycloid(A, B, n=y_smooth.size)\n",
      "\n",
      "# Generate a parabola\n",
      "y1 = (x- B[0]*0.75)**2\n",
      "\n",
      "# Generate a wavey curve\n",
      "y2 = 0.5*x + 2.5*np.sin(0.1*x + 1*np.pi/2) \n",
      "\n",
      "# Now fit everything together between A and B\n",
      "y0, y_smooth, y1, y2, y3 = fit_multiple([y0, y_smooth,y1, y2, x], A, B)\n",
      "\n",
      "t_drawn = total_time(x, y_smooth)\n",
      "t_wavey = total_time(x, y2)\n",
      "t_parabola = total_time(x, y1)\n",
      "t_straight = total_time(x, y3)\n",
      "t_cycloid = total_time(x, y0)\n",
      "\n",
      "t_max = max([t_drawn, t_wavey, t_parabola, t_straight, t_cycloid])\n",
      "\n",
      "# plots\n",
      "names = [\"Drawn\",\"Wavey\",\"Parabola\",\"Straight\",\"Inv. Cycloid\"]\n",
      "Ys = [y_smooth, y2, y1, y3, y0]\n",
      "\n",
      "def animate(t):\n",
      "    Animate_frame(Ys, x, t, names)\n",
      "    text(A[0]+5, A[1]+10, \"A\", fontsize=20)\n",
      "    text(B[0]+5, B[1]+10, \"B\", fontsize=20)\n",
      "\n",
      "t_slider = FloatSliderWidget(min=0, max=t_max+0.001, step=t_max/100, value=0)\n",
      "w=interactive(animate, t=t_slider)\n",
      "display(w)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print\n",
      "print \"Estimated total descent times:\"\n",
      "print 30*\"-\"\n",
      "print \"Drawn :   \", t_drawn\n",
      "print \"Wavey :   \", t_wavey\n",
      "print \"Parabola :\", t_parabola\n",
      "print \"Straight :\", t_straight\n",
      "print \"Cycloid : \", t_cycloid\n",
      "print"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Finally, generate a 3D STL file of marble ramps based on the above curves"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# -- this will generate the raw OpenSCAD code that will create the marble ramps\n",
      "#    unless re-scaling in the 3D slicing software before printing, all units are in mm\n",
      "\n",
      "name = \"my_ramps\" # name to use for the OpenSCAD modules and stl file\n",
      "\n",
      "thickness = 8        # thickness of the base\n",
      "length = B[0] - A[0] # length of the ramps\n",
      "n = 250              # resolution of the ramp surfaces\n",
      "back_height = 10      # height of the backing on the starting (A) side\n",
      "back_depth = 15      # depth of the backing on the starting (A) side\n",
      "channel_width = 20   # width of the marble channels\n",
      "channel_depth = 10    # depth of the marble channels\n",
      "wall_width = 1.5       # thickness of the channel walls\n",
      "endstop_thickness = 5 # thickness of the end stop block\n",
      "endstop_height = 30  # height of the end stop block\n",
      "\n",
      "code = ramp([y_smooth,y2, y1, y3, y0], \n",
      "            thickness=thickness, \n",
      "            length=length, \n",
      "            start_backdepth = back_depth,\n",
      "            endstop_thickness=endstop_thickness,\n",
      "            endstop_height=endstop_height,\n",
      "            channel_width = channel_width,\n",
      "            wall_width = wall_width,\n",
      "            channel_depth = channel_depth,\n",
      "            n=n, \n",
      "            start_backheight=back_height, \n",
      "            name=\"ramp\")\n",
      "\n",
      "# Now, save the scad file, and export to .stl from the command line (may take a minute)\n",
      "render_stl(openscad_path, code, name)\n",
      "\n",
      "# Show the STL file below \n",
      "show_stl(name)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The last steps above (STL generation and preview) may take a minute or so to complete.\n",
      "\n",
      "## Notes\n",
      "\n",
      "The ramps seen are, starting on the left (facing the object from the front end/lower part):\n",
      "\n",
      "- Your drawn curve\n",
      "- Wavey curve\n",
      "- Parabola\n",
      "- Straight line\n",
      "- Inverted cycloid\n",
      "\n",
      "If you print out these ramps and test them with marbles, you'll find that a marble rolling down the last curve (the inverted cycloid) will tend to reach the end point the fastest. \n",
      "\n",
      "When prepping the STL file for printing, you can scale the object to fit the volume of your printer **as long as you scale each axis uniformly (together)**. If the length vs height ratio of the ramps is modified, then the descent properties of the curve may be very different than what was computed in this notebook.\n",
      "\n",
      "## Derivation of the inverted cycloid as the optimal solution\n",
      "Neglecting friction and drag effects, the total time for an object to roll from A$=\\left(x_0, y_0\\right)$ to B$=\\left(x_1, y_1\\right)$ along a curve $y\\left(x\\right)$ in a uniform graviational field can be computed as\n",
      "\n",
      "$$T = \\int_{x_0}^{x_1} \\left( \\left(y^\\prime\\left(x\\right) \\right)^2 + 1 \\right)^{\\frac{1}{2}} \\left(2 g y\\left(x\\right) \\right)^{-\\frac{1}{2}} dx$$\n",
      "\n",
      "where $y'\\left(x\\right) = \\frac{dy}{dx}$ and $g$ is a gravitational acceleration constant. For simplicity, let $x_0 = y_0 = 0$ and let $x_1$ and $y_1$ shift appropriately to $x_n,y_n$ in order to represent the same total length and change of height.\n",
      "So the Brachistochrone problem then has the form:\n",
      "\n",
      "\n",
      "\n",
      "$$\\min_{y\\left(x\\right) \\in C^2[0,x_1]}  \\int_{0}^{x_1} \\left( \\left(y' \\right)^2 + 1 \\right)^{\\frac{1}{2}} \\left(2 g y \\right)^{-\\frac{1}{2}} dx$$\n",
      "\n",
      "$$ =\\sqrt{2g }^{-1} \\min_{y\\left(x\\right) \\in C^2[0,x_1]} \\int_{0}^{x_1} F\\left(y'\\left(x\\right), y\\left(x\\right), x\\right)  dx $$\n",
      "\n",
      "subject to the boundary conditions $y\\left(0 \\right ) = 0$ and $y\\left(x_n \\right ) = y_n$, where\n",
      "$$F\\left(y'\\left(x\\right), y\\left(x\\right), x\\right) = \\sqrt{\\frac{ y'\\left(x \\right)^2 + 1 }{y}}$$\n",
      "\n",
      "[The Calculus of Variations](http://en.wikipedia.org/wiki/Calculus_of_variations) is a subject concerned with solving optimization problems of this type, where the values of an entire function between two points collectively serve as the parameter to be varied.\n",
      "\n",
      "As it turns out, a first-order necessary condition for a local minimizer $y$ is that it satisifies the associated [Euler-Lagrange equation](http://en.wikipedia.org/wiki/Euler%E2%80%93Lagrange_equation):\n",
      "\n",
      "$$\\frac{\\partial F}{\\partial y} - \\frac{d}{dx} \\frac{\\partial F}{\\partial y'} = 0$$\n",
      "\n",
      "and if $F$ has no explicit dependance on $x$, [then the E-L equation takes the form](http://en.wikipedia.org/wiki/Beltrami_identity)\n",
      "\n",
      "$$F - y'\\frac{\\partial F}{\\partial y'} = C  $$\n",
      "\n",
      "for a constant $C$. For the Brachistochrone problem\n",
      "\n",
      "$$ y'\\frac{\\partial F}{\\partial y'} = y'^2 \\left(y\\left(1+y'^2\\right) \\right)^{-\\frac{1}{2}} $$\n",
      "\n",
      "which leads to\n",
      "\n",
      "$$ \\sqrt{\\frac{ y'^2 + 1 }{y}} -  y'^2 \\left(y\\left(1+y'^2\\right) \\right)^{-\\frac{1}{2}} = C $$\n",
      "\n",
      "$$\\Rightarrow \\frac{1}{\\sqrt{y \\left(1 + y'^2 \\right)}} = C$$\n",
      "\n",
      "$$\\Rightarrow y \\left(1 + y'^2 \\right) = \\frac{1}{C^2}$$\n",
      "\n",
      "which is solved parametrically by the equations\n",
      "\n",
      "$$x = \\frac{1}{2C^2}\\left(\\theta - \\sin\\theta \\right) = a\\left(\\theta - \\sin\\theta \\right) $$\n",
      "$$y = \\frac{1}{2C^2}\\left(1 - \\cos\\theta \\right) =  a\\left(1 - \\cos\\theta \\right) $$\n",
      "\n",
      "for a constant $a$, [which are recognized as the equations of a cycloid](http://en.wikipedia.org/wiki/Cycloid#Equations).\n",
      "\n",
      "Shifting the problem back to A$=\\left(x_0, y_0\\right)$ to B$=\\left(x_1, y_1\\right)$, the equations are\n",
      "\n",
      "$$x = x_0 + a\\left(\\theta - \\sin\\theta \\right) $$\n",
      "$$y = y_0 + a\\left(1 - \\cos\\theta \\right) $$\n",
      "\n",
      "The constant $a$ can be computed based on known information. First, we know that at the end point B\n",
      "\n",
      "$$\\frac{y_0 - y_1}{x_1 - x_0} =\\frac{\\left(1 - \\cos\\theta \\right)}{\\left(\\theta - \\sin\\theta \\right)}$$\n",
      "\n",
      "which we can solve this numerically for $\\theta \\in [0, 2\\pi]$. Then we can plugin $\\theta$ into \n",
      "\n",
      "$$x = x_0 + a\\left(\\theta - \\sin\\theta \\right) $$\n",
      "\n",
      "to determine $a$.\n",
      "\n",
      "The `cycloid` function in the file `lib/cycloid.py` computes a Cycloid between the two points given in the first cell of this notebook in exactly this manner. "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}