Python for the Dynamicist

85bba1ca66eb909a289448a90e88f53a?s=47 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

85bba1ca66eb909a289448a90e88f53a?s=128

Jim Crist

March 11, 2015
Tweet

Transcript

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

  2. Overview • A brief Python Overview • Introduction to SymPy

    • Multibody Dynamics Review • Examples!
  3. What is Python? “A multiple paradigm, dynamically typed, interpreted, high

    level language focusing on readability.”
  4. Core Data Types • Integers 1, 2, 3 • Floats

    1.1, 2.5, 3.14 • Strings "apple", 'apple' • Constants True, False, None
  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}
  6. Functions def count_chars(string): """Return the number of characters in a

    string.""" count = 0 for char in string: count += 1 return count
  7. Why use Python??? • Free • General purpose language •

    Great Standard Library • Fast enough (usually) • Easily extensible • Active community
  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!
  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
  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]])
  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
  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]])
  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
  14. </python evangelism>

  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
  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
  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)
  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)]
  19. Submodules (just a few…) • Calculus • Linear Algebra •

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

    Set Theory • Combinatorics • Differential Geometry • Tensors • Classical Mechanics • etc… What this talk is about
  21. A VERY BRIEF OVERVIEW OF CLASSICAL MECHANICS

  22. Classical Mechanics • Describes the motion of macroscopic objects •

    Concerned with solving = = −1 ∙ • Many ways to do so: • Newton • Lagrange • Kane
  23. Vectors & Reference Frames

  24. Generalized Coordinates & Speeds

  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)
  26. Constraints • Pendulum has 1 DOF, but 2 coordinates •

    Length is always • 1 , 2 , = 1 2 + 2 2 − 2 = 0
  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
  28. Equations of Motion ℎ , = 0×1 ℎ , ,

    = 0×1 , , , = 0×1 , , , , , = 0 −+ ×1 with , ∈ ℝ , ∈ ℝ ∈ ℝ ∈ ℝ
  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
  30. Equations of Motion Can be rearranged to form: (, )

    = (, , , ) or = , −1(, , , )
  31. Classical Mechanics • Inverted Pendulum, 1 Link

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

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

    2.5 font
  34. DERIVE EQUATIONS IN SYMPY Solution:

  35. github.com/pydy/pydy

  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
  37. Demo (try it at github.com/jcrist/talks/tree/master/dynamics_talk)

  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 ()
  39. Example – Single Link Pendulum • Derived EOM for minimal,

    and nonminimal formulations
  40. Example – Single Link Pendulum Response for small angle deviation

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

    - simple jacobian linearization Wrong…
  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)
  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
  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)
  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
  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
  47. Questions???