code/chap25.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling and Simulation in Python\n",
    "\n",
    "Chapter 25\n",
    "\n",
    "Copyright 2017 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
    "\n",
    "# import functions from the modsim.py module\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Teapots and Turntables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tables in Chinese restaurants often have a rotating tray or turntable that makes it easy for customers to share dishes.  These turntables are supported by low-friction bearings that allow them to turn easily and glide.  However, they can be heavy, especially when they are loaded with food, so they have a high moment of inertia.\n",
    "\n",
    "Suppose I am sitting at a table with a pot of tea on the turntable directly in front of me, and the person sitting directly opposite asks me to pass the tea.  I push on the edge of the turntable with 1 Newton of force until it has turned 0.5 radians, then let go.  The turntable glides until it comes to a stop 1.5 radians from the starting position.  How much force should I apply for a second push so the teapot glides to a stop directly opposite me?\n",
    "\n",
    "The following figure shows the scenario, where `F` is the force I apply to the turntable at the perimeter, perpendicular to the moment arm, `r`, and `tau` is the resulting torque.  The blue circle near the bottom is the teapot.\n",
    "\n",
    "![](diagrams/teapot.png)\n",
    "\n",
    "We'll answer this question in these steps:\n",
    "\n",
    "1.  We'll use the results from the first push to estimate the coefficient of friction for the turntable.\n",
    "\n",
    "2.  We'll use that coefficient of friction to estimate the force needed to rotate the turntable through the remaining angle.\n",
    "\n",
    "Our simulation will use the following parameters:\n",
    "\n",
    "1.  The radius of the turntable is 0.5 meters, and its weight is 7 kg.\n",
    "\n",
    "2.  The teapot weights 0.3 kg, and it sits 0.4 meters from the center of the turntable.\n",
    "\n",
    "As usual, I'll get units from Pint."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "radian = UNITS.radian\n",
    "m = UNITS.meter\n",
    "s = UNITS.second\n",
    "kg = UNITS.kilogram\n",
    "N = UNITS.newton"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And store the parameters in a `Params` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "params = Params(radius_disk=0.5*m,\n",
    "                mass_disk=7*kg,\n",
    "                radius_pot=0.4*m,\n",
    "                mass_pot=0.3*kg,\n",
    "                force=1*N,\n",
    "                torque_friction=0.2*N*m,\n",
    "                theta_end=0.5*radian)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`make_system` creates the initial state, `init`, and computes the total moment of inertia for the turntable and the teapot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def make_system(params):\n",
    "    \"\"\"Make a system object.\n",
    "    \n",
    "    params: Params object\n",
    "    \n",
    "    returns: System object\n",
    "    \"\"\"\n",
    "    unpack(params)\n",
    "    \n",
    "    init = State(theta=0, omega=0)\n",
    "    \n",
    "    I_disk = mass_disk * radius_disk**2 / 2\n",
    "    I_pot = mass_pot * radius_pot**2\n",
    "    \n",
    "    return System(params, init=init, t_end=20*s,\n",
    "                  I=I_disk+I_pot)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the `System` object we'll use for the first phase of the simulation, while I am pushing the turntable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "system1 = make_system(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulation\n",
    "\n",
    "When I stop pushing on the turntable, the angular acceleration changes abruptly.  We could implement the slope function with an `if` statement that checks the value of `theta` and sets `force` accordingly.  And for a coarse model like this one, that might be fine.  But we will get more accurate results if we simulate the system in two phases:\n",
    "\n",
    "1.  During the first phase, force is constant, and we run until `theta` is 0.5 radians.\n",
    "\n",
    "2.  During the second phase, force is 0, and we run until `omega` is 0.\n",
    "\n",
    "Then we can combine the results of the two phases into a single `TimeFrame`.\n",
    "\n",
    "Here's the slope function we'll use:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def slope_func(state, t, system):\n",
    "    \"\"\"Computes the derivatives of the state variables.\n",
    "    \n",
    "    state: State object\n",
    "    t: time\n",
    "    system: System object \n",
    "    \n",
    "    returns: sequence of derivatives\n",
    "    \"\"\"\n",
    "    theta, omega = state\n",
    "    unpack(system)\n",
    "    \n",
    "    torque = radius_disk * force - torque_friction\n",
    "    alpha = torque / I\n",
    "    \n",
    "    return omega, alpha "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As always, we'll test the slope function before running the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope_func(system1.init, 0*s, system1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's an event function that stops the simulation when `theta` reaches `theta_end`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def event_func1(state, t, system):\n",
    "    \"\"\"Stops when theta reaches theta_end.\n",
    "    \n",
    "    state: State object\n",
    "    t: time\n",
    "    system: System object \n",
    "    \n",
    "    returns: difference from target\n",
    "    \"\"\"\n",
    "    theta, omega = state\n",
    "    unpack(system)\n",
    "    return theta - theta_end "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can run the first phase."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "results1, details1 = run_ode_solver(system1, slope_func,\n",
    "                                    events=event_func1, max_step=0.1*s)\n",
    "details1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And look at the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "results1.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Phase 2\n",
    "\n",
    "Before we run the second phase, we have to extract the final time and state of the first phase."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_0 = get_last_label(results1) * s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And make an initial `State` object for Phase 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta, omega = get_last_value(results1)\n",
    "init2 = State(theta=theta*radian, omega=omega*radian/s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And a new `System` object with `force=0`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "system2 = System(system1, t_0=t_0, init=init2, force=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's an event function that stops when angular velocity is 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def event_func2(state, t, system):\n",
    "    \"\"\"Stops when omega is 0.\n",
    "    \n",
    "    state: State object\n",
    "    t: time\n",
    "    system: System object \n",
    "    \n",
    "    returns: omega\n",
    "    \"\"\"\n",
    "    theta, omega = state\n",
    "    return omega"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can run the second phase."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "results2, details2 = run_ode_solver(system2, slope_func,\n",
    "                         events=event_func2, max_step=0.1*s)\n",
    "details2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And check the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "results2.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pandas provides `combine_first`, which combines `results1` and `results2`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = results1.combine_first(results2)\n",
    "results.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can plot `theta` for both phases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_theta(results):\n",
    "    plot(results.theta, label='theta')\n",
    "    decorate(xlabel='Time (s)',\n",
    "             ylabel='Angle (rad)')\n",
    "    \n",
    "plot_theta(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And `omega`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_omega(results):\n",
    "    plot(results.omega, label='omega', color='C1')\n",
    "    decorate(xlabel='Time (s)',\n",
    "             ylabel='Angular velocity (rad/s)')\n",
    "    \n",
    "plot_omega(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "subplot(2, 1, 1)\n",
    "plot_theta(results)\n",
    "subplot(2, 1, 2)\n",
    "plot_omega(results)\n",
    "savefig('figs/chap25-fig01.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimating friction\n",
    "\n",
    "Let's take the code from the previous section and wrap it in a function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_two_phases(force, torque_friction, params):\n",
    "    \"\"\"Run both phases.\n",
    "    \n",
    "    force: force applied to the turntable\n",
    "    torque_friction: friction due to torque\n",
    "    params: Params object\n",
    "    \n",
    "    returns: TimeFrame of simulation results\n",
    "    \"\"\"\n",
    "    # put the specified parameters into the Params object\n",
    "    params = Params(params, force=force, torque_friction=torque_friction)\n",
    "\n",
    "    # run phase 1\n",
    "    system1 = make_system(params)\n",
    "    results1, details1 = run_ode_solver(system1, slope_func, \n",
    "                             events=event_func1, max_step=0.1*s)\n",
    "\n",
    "    # get the final state from phase 1\n",
    "    t_0 = get_last_label(results1) * s\n",
    "    theta, omega = get_last_value(results1)\n",
    "    init2 = State(theta=theta, omega=omega)\n",
    "    \n",
    "    # run phase 2\n",
    "    system2 = System(system1, t_0=t_0, init=init2, force=0)\n",
    "    results2, details2 = run_ode_solver(system2, slope_func, \n",
    "                             events=event_func2, max_step=0.1*s)\n",
    "    \n",
    "    # combine and return the results\n",
    "    results = results1.combine_first(results2)\n",
    "    return results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's test it with the same parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "force = 1*N\n",
    "torque_friction = 0.2*N*m\n",
    "results = run_two_phases(force, torque_friction, params)\n",
    "results.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And check the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta_final = get_last_value(results.theta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To estimate the coefficient of friction, we'll use `fsolve`, which doesn't always work well with units.\n",
    "\n",
    "So for the rest of this example, we'll run without units.  Here's a version of the `Params` object with no units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "params_nodim = remove_units(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the error function we'll use with `fsolve`.  It takes a hypothetical value for `torque_friction` and returns the difference between `theta_final` and the observed result of the first push."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "def error_func1(torque_friction, params):\n",
    "    \"\"\"Error function for fsolve.\n",
    "    \n",
    "    torque_friction: hypothetical value\n",
    "    params: Params object\n",
    "    \n",
    "    returns: offset from target value\n",
    "    \"\"\"\n",
    "    force = 1\n",
    "    results = run_two_phases(force, torque_friction, params)\n",
    "    theta_final = get_last_value(results.theta)\n",
    "    print(torque_friction, theta_final)\n",
    "    return theta_final - 1.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Testing the error function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "guess = 0.2\n",
    "error_func1(guess, params_nodim)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And running `fsolve`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = fsolve(error_func1, guess, params_nodim)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is an array, so we'll extract the first (and only) element."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "torque_friction = res[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And do a test run with the estimated value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "force = 1\n",
    "results = run_two_phases(force, torque_friction, params_nodim)\n",
    "theta_final = get_last_value(results.theta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks good.\n",
    "\n",
    "### Exercises\n",
    "\n",
    "Now finish off the example by estimating the force that delivers the teapot to the desired position.\n",
    "\n",
    "Write an error function that takes `force` and `params` and returns the offset from the desired angle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test the error function with `force=1`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And run `fsolve` to find the desired force."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "force = res[0]\n",
    "results = run_two_phases(force, torque_friction, params_nodim)\n",
    "theta_final = get_last_value(results.theta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Now suppose my friend pours 0.1 kg of tea and puts the teapot back on the turntable at distance 0.3 meters from the center.  If I ask for the tea back, how much force should they apply, over an arc of 0.5 radians, to make the teapot glide to a stop back in front of me?  You can assume that torque due to friction is proportional to the total mass of the teapot and the turntable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "subplot(2, 1, 1)\n",
    "plot_theta(results)\n",
    "subplot(2, 1, 2)\n",
    "plot_omega(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "**Exercise:** Run these simulations with different values for the weight and size of the turntable, or weight and location of the teapot.  What effect do they have on the results?\n",
    "\n",
    "If you do some analysis, you might find that you can answer the original question without doing any simulation."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}