Multibody Dynamics with SymPy and PyDy

A talk I gave for my research lab about PyDy [1], SymPy [2], and my experience in Google Summer of Code [3].

Demo Notebook available on github: https://github.com/jcrist/talks/tree/master/pydy_talk.

[1] https://github.com/pydy/pydy
[2] https://github.com/sympy/sympy
[3] https://www.google-melange.com/gsoc/homepage/google/gsoc2014

Transcript

1. Multibody Dynamics with
   SymPy and PyDy
   Jim Crist
   (Or what I did this summer in 30 slides or less)
   

   View Slide

2. 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

3. 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

4. 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

5. 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

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

   View Slide

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

   View Slide

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

   View Slide

9. A VERY BRIEF OVERVIEW OF
   MULTIBODY
   DYNAMICS
   

   View Slide

10. Vectors & Reference Frames
    
    
    
    
    
    
    
    
    
    

    View Slide

11. Vectors & Reference Frames
    

    View Slide

12. Generalized Coordinates & Speeds
    

    View Slide

13. 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

14. Constraints
    • Pendulum has 1 DOF, but 2
    coordinates
    • Length is always
    •
    1
    , 2
    , = 1
    2 + 2
    2 − 2 = 0
    

    View Slide

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

    View Slide

16. Equations of Motion
    
    , = 0×1
    
    , , = 0×1
    
    , , , = 0×1
    
    , , , = 0×1
    
    , , , ,, = 0 0−+ ×1
    with
    , ∈ ℝ
    , ∈ ℝ
    ∈ ℝ
    ∈ ℝ
    

    View Slide

17. Equations of Motion
    Can be rearranged to form:
    (, )
    
    
    
    = (, , , ,)
    or
    
    
    
    = , −1(, , ,, )
    

    View Slide

18. Classical Mechanics
    • Inverted Pendulum, 1 Link
    

    View Slide

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

    View Slide

20. DERIVE
    EQUATIONS
    IN SYMPY
    Solution:
    

    View Slide

21. github.com/pydy/pydy
    

    View Slide

22. Stuff I Did…
    

    View Slide

23. 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

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

    View Slide

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

    View Slide

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

    View Slide

27. Linearization
    • Derived and implemented algorithm 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

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

    View Slide

29. Other Stuff
    • Rewrote solver – 100x speed boost for (some) benchmarks
    • Improvements to substitution routines
    • Now () not (2)
    • Bug fixes for lots of things
    • Documentation
    

    View Slide

30. Demo
    (try it at github.com/jcrist/talks/tree/master/pydy_talk)
    

    View Slide

31. 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
    

    View Slide

32. Future Work
    • O(n) articulated methods
    • Featherstone/Jain
    • Finish work on LAPACK/BLAS/ODEPACK integration
    • Will remove python layer bottleneck for simulation
    • Integration with CasADi
    

    View Slide