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

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