Skip to content

Commit

Permalink
Merge pull request #146 from tkoolen/sympy
Browse files Browse the repository at this point in the history
Add symbolic dynamics example using SymPy.jl
  • Loading branch information
tkoolen authored Jan 31, 2017
2 parents 32cf3bb + 83a448b commit f2bf264
Show file tree
Hide file tree
Showing 12 changed files with 373 additions and 42 deletions.
3 changes: 2 additions & 1 deletion .travis.yml
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@ notifications:
before_install:
- if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
script:
- julia -e 'Pkg.clone(pwd()); Pkg.add("Rotations"); Pkg.checkout("Rotations"); Pkg.build("RigidBodyDynamics"); Pkg.test("RigidBodyDynamics"; coverage=true)'
- julia -e 'ENV["PYTHON"]=""; Pkg.clone(pwd()); Pkg.build("RigidBodyDynamics"); Pkg.test("RigidBodyDynamics"; coverage=true)'
# Note: PYTHON env is to ensure that PyCall uses the Conda.jl package's Miniconda distribution within Julia. Otherwise the sympy Python module won't be installed/imported properly.
after_success:
- julia -e 'cd(Pkg.dir("RigidBodyDynamics")); Pkg.add("Coverage"); using Coverage; Codecov.submit(Codecov.process_folder())'
329 changes: 329 additions & 0 deletions examples/Symbolic double pendulum.ipynb
Original file line number Diff line number Diff line change
@@ -0,0 +1,329 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"using RigidBodyDynamics\n",
"using StaticArrays\n",
"using SymPy"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Create symbolic parameters\n",
"* Masses: $m_1, m_2$\n",
"* Mass moments of inertia (about center of mass): $I_1, I_2$\n",
"* Link lengths: $l_1, l_2$\n",
"* Center of mass locations (w.r.t. preceding joint axis): $c_1, c_2$\n",
"* Gravitational acceleration: $g$"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/latex": [
"\\begin{bmatrix}m_{1}\\\\m_{2}\\\\I_{1}\\\\I_{2}\\\\l_{1}\\\\l_{2}\\\\c_{1}\\\\c_{2}\\\\g\\end{bmatrix}"
],
"text/plain": [
"9-element Array{SymPy.Sym,1}\n",
"[m_1]\n",
"[ ]\n",
"[m_2]\n",
"[ ]\n",
"[I_1]\n",
"[ ]\n",
"[I_2]\n",
"[ ]\n",
"[l_1]\n",
"[ ]\n",
"[l_2]\n",
"[ ]\n",
"[c_1]\n",
"[ ]\n",
"[c_2]\n",
"[ ]\n",
"[ g ]"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"inertias = @syms m_1 m_2 I_1 I_2 positive = true\n",
"lengths = @syms l_1 l_2 c_1 c_2 real = true\n",
"gravitationalAcceleration = @syms g real = true\n",
"params = [inertias..., lengths..., gravitationalAcceleration...]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Create double pendulum `Mechanism`\n",
"\n",
"A `Mechanism` contains the joint layout and inertia parameters, but no state information."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"Vertex: world (root)\n",
" Vertex: upper_link, Edge: shoulder\n",
" Vertex: lower_link, Edge: elbow"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"T = Sym # the 'scalar type' of the Mechanism we'll construct\n",
"axis = SVector(zero(T), one(T), zero(T)) # axis of rotation for each of the joints\n",
"doublePendulum = Mechanism(RigidBody{T}(\"world\"); gravity = SVector(0, 0, g))\n",
"world = root_body(doublePendulum) # the fixed 'world' rigid body\n",
"\n",
"# Attach the first (upper) link to the world via a revolute joint named 'shoulder'\n",
"inertia1 = SpatialInertia(CartesianFrame3D(\"upper_link\"), I_1 * axis * axis', SVector(0, 0, c_1), m_1)\n",
"body1 = RigidBody(inertia1)\n",
"joint1 = Joint(\"shoulder\", Revolute(axis))\n",
"joint1ToWorld = Transform3D{T}(joint1.frameBefore, default_frame(world))\n",
"attach!(doublePendulum, world, joint1, joint1ToWorld, body1)\n",
"\n",
"# Attach the second (lower) link to the world via a revolute joint named 'elbow'\n",
"inertia2 = SpatialInertia(CartesianFrame3D(\"lower_link\"), I_2 * axis * axis', SVector(0, 0, c_2), m_2)\n",
"body2 = RigidBody(inertia2)\n",
"joint2 = Joint(\"elbow\", Revolute(axis))\n",
"joint2ToBody1 = Transform3D(joint2.frameBefore, default_frame(body1), SVector(0, 0, l_1))\n",
"attach!(doublePendulum, body1, joint2, joint2ToBody1, body2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Create `MechanismState` associated with the double pendulum `Mechanism`\n",
"\n",
"A `MechanismState` stores all state-dependent information associated with a `Mechanism`."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x = MechanismState(T, doublePendulum);"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/latex": [
"\\begin{bmatrix}q_{1}\\\\q_{2}\\end{bmatrix}"
],
"text/plain": [
"2-element Array{SymPy.Sym,1}\n",
"[q_1]\n",
"[ ]\n",
"[q_2]"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Set the joint configuration vector of the MechanismState to a new vector of symbolic variables\n",
"q = configuration_vector(x)[:] = [symbols(\"q_$i\", real = true) for i = 1 : num_positions(x)]"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/latex": [
"\\begin{bmatrix}v_{1}\\\\v_{2}\\end{bmatrix}"
],
"text/plain": [
"2-element Array{SymPy.Sym,1}\n",
"[v_1]\n",
"[ ]\n",
"[v_2]"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Set the joint velocity vector of the MechanismState to a new vector of symbolic variables\n",
"v = velocity_vector(x)[:] = [symbols(\"v_$i\", real = true) for i = 1 : num_positions(x)]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Compute dynamical quantities in symbolic form"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/latex": [
"\\begin{bmatrix}I_{1} + I_{2} + 2 c_{2} l_{1} \\cos{\\left (q_{2} \\right )} + l_{1}^{2} m_{2}&I_{2} + c_{2} l_{1} \\cos{\\left (q_{2} \\right )}\\\\I_{2} + c_{2} l_{1} \\cos{\\left (q_{2} \\right )}&I_{2}\\end{bmatrix}"
],
"text/plain": [
"2×2 Array{SymPy.Sym,2}\n",
"[ 2 ]\n",
"[I_1 + I_2 + 2*c_2*l_1*cos(q_2) + l_1 *m_2 I_2 + c_2*l_1*cos(q_2)]\n",
"[ ]\n",
"[ I_2 + c_2*l_1*cos(q_2) I_2 ]"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Mass matrix\n",
"M = mass_matrix(x)\n",
"function simplify!(a::Array{SymPy.Sym})\n",
" for i in eachindex(a)\n",
" a[i] = simplify(a[i])\n",
" end\n",
" a\n",
"end\n",
"simplify!(M.data) # Note: M is a Symmetric matrix type; need to simplify the underlying data\n",
"full(M) # convert to full form so that it is pretty-printed (minor bug in SymPy.jl)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/latex": [
"$$\\frac{I_{1} v_{1}^{2}}{2} + \\frac{I_{2} v_{1}^{2}}{2} + I_{2} v_{1} v_{2} + \\frac{I_{2} v_{2}^{2}}{2} + c_{2} l_{1} v_{1}^{2} \\cos{\\left (q_{2} \\right )} + c_{2} l_{1} v_{1} v_{2} \\cos{\\left (q_{2} \\right )} + \\frac{l_{1}^{2} m_{2}}{2} v_{1}^{2}$$"
],
"text/plain": [
" 2 2 2 \n",
"I_1*v_1 I_2*v_1 I_2*v_2 2 \n",
"-------- + -------- + I_2*v_1*v_2 + -------- + c_2*l_1*v_1 *cos(q_2) + c_2*l_1\n",
" 2 2 2 \n",
"\n",
" 2 2\n",
" l_1 *m_2*v_1 \n",
"*v_1*v_2*cos(q_2) + -------------\n",
" 2 "
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Kinetic energy\n",
"simplify(kinetic_energy(x))"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/latex": [
"$$- g \\left(c_{1} \\cos{\\left (q_{1} \\right )} + c_{2} \\cos{\\left (q_{1} + q_{2} \\right )} + l_{1} m_{2} \\cos{\\left (q_{1} \\right )}\\right)$$"
],
"text/plain": [
"-g*(c_1*cos(q_1) + c_2*cos(q_1 + q_2) + l_1*m_2*cos(q_1))"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Potential energy\n",
"simplify(potential_energy(x))"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 0.5.0",
"language": "julia",
"name": "julia-0.5"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "0.5.0"
}
},
"nbformat": 4,
"nbformat_minor": 1
}
10 changes: 5 additions & 5 deletions src/frames.jl
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@ end


# Transform between frames
immutable Transform3D{T<:Real}
immutable Transform3D{T<:Number}
from::CartesianFrame3D
to::CartesianFrame3D
rot::RotMatrix3{T}
Expand Down Expand Up @@ -105,14 +105,14 @@ for VectorType in (:FreeVector3D, :Point3D)
end

$VectorType{V}(frame::CartesianFrame3D, v::V) = $VectorType{V}(frame, v)
$VectorType{T<:Real}(::Type{T}, frame::CartesianFrame3D) = $VectorType(frame, zeros(SVector{3, T}))
$VectorType{T<:Number}(::Type{T}, frame::CartesianFrame3D) = $VectorType(frame, zeros(SVector{3, T}))

convert{V}(::Type{$VectorType{V}}, p::$VectorType{V}) = p
convert{V}(::Type{$VectorType{V}}, p::$VectorType) = $VectorType(p.frame, convert(V, p.v))

(/){S<:Real}(p::$VectorType, s::S) = $VectorType(p.frame, p.v / s)
(*){S<:Real}(p::$VectorType, s::S) = $VectorType(p.frame, p.v * s)
(*){S<:Real}(s::S, p::$VectorType) = $VectorType(p.frame, s * p.v)
(/){S<:Number}(p::$VectorType, s::S) = $VectorType(p.frame, p.v / s)
(*){S<:Number}(p::$VectorType, s::S) = $VectorType(p.frame, p.v * s)
(*){S<:Number}(s::S, p::$VectorType) = $VectorType(p.frame, s * p.v)

rand{T}(::Type{$VectorType}, ::Type{T}, frame::CartesianFrame3D) = $VectorType(frame, rand(SVector{3, T}))
show(io::IO, p::$VectorType) = print(io, "$($(VectorType).name.name) in \"$(name(p.frame))\": $(p.v)")
Expand Down
4 changes: 2 additions & 2 deletions src/joint.jl
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
type Joint{T<:Real}
type Joint{T<:Number}
name::String
frameBefore::CartesianFrame3D
frameAfter::CartesianFrame3D
Expand All @@ -7,7 +7,7 @@ type Joint{T<:Real}
Joint(name::String, jointType::JointType{T}) = new(name, CartesianFrame3D(string("before_", name)), CartesianFrame3D(string("after_", name)), jointType)
end

Joint{T<:Real}(name::String, jointType::JointType{T}) = Joint{T}(name, jointType)
Joint{T<:Number}(name::String, jointType::JointType{T}) = Joint{T}(name, jointType)

show(io::IO, joint::Joint) = print(io, "Joint \"$(joint.name)\": $(joint.jointType)")
showcompact(io::IO, joint::Joint) = print(io, "$(joint.name)")
Expand Down
Loading

0 comments on commit f2bf264

Please sign in to comment.