Python for the Dynamicist Jim Crist 3/11/15

Overview • A brief Python Overview • Introduction to SymPy • Multibody Dynamics Review • Examples!

What is Python? “A multiple paradigm, dynamically typed, interpreted, high level language focusing on readability.”

Core Data Types • Integers 1, 2, 3 • Floats 1.1, 2.5, 3.14 • Strings "apple", 'apple' • Constants True, False, None

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}

Functions def count_chars(string): """Return the number of characters in a string.""" count = 0 for char in string: count += 1 return count

Why use Python??? • Free • General purpose language • Great Standard Library • Fast enough (usually) • Easily extensible • Active community

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!

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

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

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

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

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

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

A VERY BRIEF OVERVIEW OF CLASSICAL MECHANICS

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

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 Lagrange Multipliers • Uses generalized forces ( and ∗) • Can be simpler to compute for large systems

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

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

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

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

DERIVE EQUATIONS IN SYMPY Solution:

github.com/pydy/pydy

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

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

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

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

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

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)

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

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

Questions???