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

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

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

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)

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)]

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

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

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

A VERY BRIEF OVERVIEW OF MULTIBODY DYNAMICS

Vectors & Reference Frames

Vectors & Reference Frames

Generalized Coordinates & Speeds

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)

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

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

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

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

Classical Mechanics • Inverted Pendulum, 1 Link

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

DERIVE EQUATIONS IN SYMPY Solution:

github.com/pydy/pydy

Stuff I Did…

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

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

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

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

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)

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

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

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

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

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