Upgrade to Pro — share decks privately, control downloads, hide ads and more …

Python for the Dynamicist

Jim Crist
March 11, 2015

Python for the Dynamicist

A presentation given for ME 8287: Advanced Robotics with Medical Applications at the University of Minnesota on using Python for rigid body dynamics analysis and control.

Materials can be found here: https://github.com/jcrist/talks/tree/master/dynamics_talk

Demo can be downloaded, or you can view it on nbviewer here: http://nbviewer.ipython.org/github/jcrist/talks/blob/master/dynamics_talk/dynamics_talk_demo.ipynb

Jim Crist

March 11, 2015
Tweet

More Decks by Jim Crist

Other Decks in Research

Transcript

  1. Python for the Dynamicist
    Jim Crist
    3/11/15

    View Slide

  2. Overview
    • A brief Python Overview
    • Introduction to SymPy
    • Multibody Dynamics Review
    • Examples!

    View Slide

  3. What is Python?
    “A multiple paradigm, dynamically typed,
    interpreted, high level language focusing
    on readability.”

    View Slide

  4. Core Data Types
    • Integers
    1, 2, 3
    • Floats
    1.1, 2.5, 3.14
    • Strings
    "apple", 'apple'
    • Constants
    True, False, None

    View Slide

  5. Core Data Types
    • Lists
    ["apple", "orange", "pear"]
    • Tuples
    ("apple", "orange", "pear")
    • Dictionaries
    {"apple": 5, "orange": 6, "pear": 4}
    • Sets
    >>> {1, 2, 3, 3}
    {1, 2, 3}

    View Slide

  6. Functions
    def count_chars(string):
    """Return the number of characters in a string."""
    count = 0
    for char in string:
    count += 1
    return count

    View Slide

  7. Why use Python???
    • Free
    • General purpose language
    • Great Standard Library
    • Fast enough (usually)
    • Easily extensible
    • Active community

    View Slide

  8. Active Community
    • Someone has probably already written the code you need
    • Python Package Index (PyPI)
    • https://pypi.python.org/pypi
    • Use pip or conda to install
    • As of 3/10/15, there were 56245 packages!

    View Slide

  9. Scientific Python Ecosystem
    • NumPy
    • Base array manipulation package
    • SciPy
    • Huge number of scientific routines
    • Wrappers for many Fortran libraries
    • Matplotlib
    • Plotting library
    • IPython
    • A better interactive environment
    • Pandas
    • Statistical analysis
    • SymPy
    • Computer algebra

    View Slide

  10. Numeric Arrays in Python
    • Not a core datatype, added with the NumPy package
    • Syntax:
    >>> from numpy import array
    >>> array([[1, 2], [3, 4]])
    array([[1, 2],
    [3, 4]])

    View Slide

  11. Numeric “Gotchas”
    • Integer division by default in Python 2
    >>> 4/3
    1
    >>> 4/3.
    1.33333333333333
    • In Python 3, automatic promotion to float
    >>> 4/3
    1.33333333333333
    >>> 4//3
    1

    View Slide

  12. Numeric “Gotchas”
    • Array multiplication is element-wise
    >>> a = array([[1, 2], [3, 4]])
    >>> a * a
    array([[1, 4],
    [9, 16]])
    • Use dot to do matrix multiplication (dot product)
    >>> a.dot(a)
    array([[7, 10],
    [15, 22]])

    View Slide

  13. To Learn More
    • Go through the Official Python Tutorial
    • Installation is easy on Mac/Linux
    • Install with system package manager (brew, apt-get, pacman…)
    • Can be a bit of a pain on Windows
    • Recommend Anaconda distribution

    View Slide


  14. View Slide

  15. What is SymPy???
    • A Computer Algebra System (CAS) written in pure Python
    • Allows you to create expressions, and do things with them
    • SymPy knows about all the math needed

    View Slide

  16. Create Expressions
    >>> from sympy import *
    # Create Symbols
    >>> a, b, c = symbols('a, b, c')
    # Combine Symbols into expressions
    >>> expr = sin(a) + cos(b**2)*c
    >>> expr
    sin(a) + cos(b**2)*c

    View Slide

  17. Manipulate Expressions
    # subs performs substitutions on expressions
    >>> expr = expr.subs(b, a**(1/2))
    >>> expr
    sin(a) + sin(a)*c
    # simplify reduces the size of expressions
    >>> expr = expr.simplify()
    >>> expr
    sin(a)*(1 + c)

    View Slide

  18. Solve Expressions
    # Create another symbol x
    >>> x = symbols('x')
    # Quadratic equation
    >>> lhs = a*x**2 + b*x + c
    >>> expr = Eq(lhs, 0)
    >>> expr
    a*x**2 + b*x + c == 0
    # Solve returns solutions to expressions
    >>> solve(expr, x)
    [(-b + sqrt(-4*a*c + b**2))/(2*a), -(b + sqrt(-4*a*c + b**2))/(2*a)]

    View Slide

  19. Submodules (just a few…)
    • Calculus
    • Linear Algebra
    • Set Theory
    • Combinatorics
    • Differential Geometry
    • Tensors
    • Classical Mechanics
    • etc…

    View Slide

  20. Submodules (just a few…)
    • Calculus
    • Linear Algebra
    • Set Theory
    • Combinatorics
    • Differential Geometry
    • Tensors
    • Classical Mechanics
    • etc…
    What this talk
    is about

    View Slide

  21. A VERY BRIEF OVERVIEW OF
    CLASSICAL
    MECHANICS

    View Slide

  22. Classical Mechanics
    • Describes the motion of macroscopic objects
    • Concerned with solving
    = = −1 ∙
    • Many ways to do so:
    • Newton
    • Lagrange
    • Kane

    View Slide

  23. Vectors & Reference Frames

    View Slide

  24. Generalized Coordinates & Speeds

    View Slide

  25. Constraints
    • If coordinates and speeds are nonminimal, need constraints
    • Holonomic
    • Relations between coordinates
    • 1
    , 2
    , 3
    , … ,
    , = 0
    • Nonholonomic
    • “Not-holonomic” (i.e. cannot be put in form above)

    View Slide

  26. Constraints
    • Pendulum has 1 DOF, but 2
    coordinates
    • Length is always

    1
    , 2
    , = 1
    2 + 2
    2 − 2 = 0

    View Slide

  27. Formulations
    • Lagrange
    • Enforces constraints with Lagrange Multipliers ()
    • Energy based formulation
    • Kane
    • Avoids need to compute Lagrange Multipliers
    • Uses generalized forces (
    and
    ∗)
    • Can be simpler to compute for large systems

    View Slide

  28. Equations of Motion

    , = 0×1

    , , = 0×1

    , , , = 0×1

    , , , , , = 0 −+ ×1
    with
    , ∈ ℝ
    , ∈ ℝ
    ∈ ℝ
    ∈ ℝ

    View Slide

  29. Equations of Motion

    , = 0×1

    , , = 0×1

    , , , = 0×1

    , , , , , = 0 −+ ×1
    with
    , ∈ ℝ
    , ∈ ℝ
    ∈ ℝ
    ∈ ℝ
    Holonomic Constraints
    Nonholonomic Constraints
    Kinematic Equations
    Dynamic Equations
    Generalized Coordinates
    Generalized Speeds
    System Inputs
    Lagrange Multipliers

    View Slide

  30. Equations of Motion
    Can be rearranged to form:
    (, )



    = (, , , )
    or



    = , −1(, , , )

    View Slide

  31. Classical Mechanics
    • Inverted Pendulum, 1 Link

    View Slide

  32. Classical Mechanics
    • Inverted Pendulum, 3 Links: Rendered in size 2.5 font

    View Slide

  33. Classical Mechanics
    • Inverted Pendulum, 3 Links: Rendered in size 2.5 font

    View Slide

  34. DERIVE
    EQUATIONS
    IN SYMPY
    Solution:

    View Slide

  35. github.com/pydy/pydy

    View Slide

  36. The PyDy Workflow
    1. Specify system symbolically
    2. Solve for EOM using either Kane’s or Lagrange’s method
    3. Analyze resulting EOM symbolically
    4. Compile equations, simulate system numerically
    5. Visualize simulation results

    View Slide

  37. Demo
    (try it at github.com/jcrist/talks/tree/master/dynamics_talk)

    View Slide

  38. Analysis - Linearization
    • For systems with constraints, need to be careful!
    • Can’t just take Jacobian, as states are interdependent
    • i.e. 1
    , , instead of 1
    ()

    View Slide

  39. Example – Single Link Pendulum
    • Derived EOM for minimal, and nonminimal formulations

    View Slide

  40. Example – Single Link Pendulum
    Response for small angle deviation - simple jacobian linearization

    View Slide

  41. Example – Single Link Pendulum
    Response for small angle deviation - simple jacobian linearization
    Wrong…

    View Slide

  42. Linearization
    • Routine implemented to robustly linearize all cases
    • Works for formulations from Kane or Lagrange
    • Handles holonomic and nonholonomic constraints
    • Applied to Whipple bicycle model [1], matches to 14 digits the
    benchmark eigenvalues
    [1] Meijaard, J.P., Papadopoulos, J.M., Ruina, A., Schwab, A.L.: Linearized dynamics
    equations for the balance and steer of a bicycle: a benchmark and review.
    Proceedings of the Royal Society A: Mathematical, Physical and Engineering
    Sciences 463(2084), 1955–1982(2007)

    View Slide

  43. Printers
    Symbolic
    Expression
    Fast
    Numeric
    Code
    Generators Wrappers
    • C
    • Fortran
    • Matlab/Octave
    • Javascript
    • Mathematica
    • Theano
    • C
    • Fortran
    • Matlab/Octave
    • Cython
    • F2Py
    • Numpy ufuncs
    Code Generation

    View Slide

  44. Caveats…
    • C/Fortran use Floating Point numbers
    • Poorly conditioned expressions can give you unrealistic results
    • No universal way to guarantee well conditioned…
    • Decent solution:
    • Simulate at several points with multiple precision numbers (GMP)
    • Repeat computation with doubles, compare for accuracy
    • Large systems -> large expressions
    • Can lead to slow derivation
    • Most systems are fine (have done 25 link inverted pendulum easily)

    View Slide

  45. Why?
    • Open Source
    • Other symbolic solvers (AutoLev*, SD/FAST) are proprietary
    • Fast
    • Slower generation than numeric methods
    • Faster simulation
    • Correct
    • No math errors for hand-calculated EOM
    • Automatic translation to C/Fortran -> no typing errors
    *Now defunct

    View Slide

  46. Future Work
    • O(n) articulated methods
    • Featherstone/Jain
    • High-level specification
    • Specify link/joint parameters to easily build up model
    • Finish work on LAPACK/BLAS/ODEPACK integration
    • Will remove Python layer bottleneck for simulation
    • Integration with Ipopt/CasADi for optimization problems

    View Slide

  47. Questions???

    View Slide