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
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???