Jim Crist
October 03, 2014
840

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

October 03, 2014

## Transcript

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

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

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

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)

5. Solve Expressions
# Create another symbol x
>>> x = symbols('x')
>>> 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)]

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

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

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

9. A VERY BRIEF OVERVIEW OF
MULTIBODY
DYNAMICS

10. Vectors & Reference Frames

11. Vectors & Reference Frames

12. Generalized Coordinates & Speeds

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)

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

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

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

16. Equations of Motion

, = 0×1

, , = 0×1

, , , = 0×1

, , , = 0×1

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

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

= (, , , ,)
or

= , −1(, , ,, )

18. Classical Mechanics

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

20. DERIVE
EQUATIONS
IN SYMPY
Solution:

21. github.com/pydy/pydy

22. Stuff I Did…

23. Linearization
• For systems with constraints, need to be careful
• Can’t just take jacobian, as states are interdependent
• i.e. 1
()

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

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

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

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)

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

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

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

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

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