Package dynkin
A toolkit for 3D dynamics and kinematics of rigid bodies using the YPR euler angle convention.
pydynkin - A python version of dynkin
A toolkit for 3D dynamics and kinematics of rigid bodies using the YPR euler angle convention.
Note: dynkin
is also available in a C++ version, available here: https://github.com/freol35241/dynkin
General
dynkin
is a set of tools for handling the dynamics and kinematics of rigid bodies in 3 dimensions (6DOFs). Some features:
- Homogenous transformation matrices
- Chained reference frames
- Idealized rigid body implementation
The fundamentals of reference frames and the kinematic relations of these are based on Theory of Applied Robotics (Reza N. Jazar) , the idealized rigid body implementation follows the outline suggested in the lectures by Fossen.
Installation
pip install dynkin
Theory intro
Some basic notions:
- A reference
Frame
is defined indynkin
as an object whichpositions
,vectors
,velocities
,accelerations
etc can be decomposed in.dynkin
represents referenceFrame
s by (4x4) homogenous transformation matrices. AFrame
is aligned (positioned, rotated) and moves (linear and angular velocity) in relation to anotherFrame
or the inertial frame (represented byNone
). - The
pose
of aFrame
is its generalized position and thetwist
of aFrame
is its generalized velocity, both in relation to the inertial frame. - All rotations in
dynkin
are internally represented by rotation matrices but the external API, so far, deals only with Euler angles of the YPR (Yaw-Pitch-Roll) convention. - A
kinematic chain
is a single-linked list ofFrame
s, where eachFrame
holds a reference to its closest parent. Any number ofFrame
s may be attached into such a chain and the chain may also have any number of branches, it is however the user´s responsibility to ensure no kinematic loops occur. - A
transform()
is an object relating twoFrame
s enabling transformation ofpositions
,vectors
,velocities
etc from oneFrame
to the other. TheFrame
s do not need to be part of the samekinematic chain
. - A
RigidBody
is a 3D body with arbitrary extent that may be described by a generalized inertia matrix (6x6). It accelerates when subject to generalized external forces (wrenches
) and rotational velocities give rise to inertial forces (coriolis and centripetal contributions).
Examples
Single frame
from dynkin import Frame, transform
frame1 = Frame(position=[1, 2, 3], attitude=[0, 0, 90], degrees=True)
# Find transformation from the inertial frame to frame1
ti1 = transform(None, frame1)
# Transformation of vector
v1_decomposed_in_frame1 = ti1.apply_vector(v1_decomposed_in_inertial_frame)
# Transformation of position
p1_decomposed_in_frame1 = ti1.apply_position(p1_decomposed_in_inertial_frame)
# Transformation of wrench
w1_decomposed_in_frame1 = ti1.apply_wrench(w1_decomposed_in_inertial_frame)
# Find the inverse transformation
t1i = ti1.inv()
# Pose of this frame, decomposed in inertial frame
frame1.get_pose()
# Twist of this frame, decomposed in inertial frame
frame.get_twist()
Two frames
from dynkin import Frame, transform
frame1 = Frame(position=[1, 2, 3], attitude=[0, 0, 90], degrees=True)
frame2 = Frame(position=[3, 2, 1], attitude=[0, 0, -90], degrees=True)
# Find transformation from frame1 to frame2
t12 = transform(frame1, frame2)
# Transformation of vector
v1_decomposed_in_frame2 = t12.apply_vector(v1_decomposed_in_frame1)
# Transformation of position
p1_decomposed_in_frame2 = t12.apply_position(p1_decomposed_in_frame1)
# Transformation of wrench
w1_decomposed_in_frame2 = t12.apply_wrench(w1_decomposed_in_frame1)
# Find the inverse transformation
t21 = t12.inv()
Kinematic chains
from dynkin import Frame, transform
frame1 = Frame(position=[1, 2, 3], attitude=[0, 0, 90], degrees=True)
frame2 = frame1.align_child(position=[3, 2, 1], attitude=[0, 0, -90], degrees=True)
frame3 = frame2.align_child(position=[1, 1, 1], attitude=[0, 0, 0], degrees=True)
# Find transformation from inertial frame to frame3
ti3 = transform(None, frame3)
# Transformation from frame3 and frame1
t31 = transform(frame3, frame1)
...
TODO: RigidBody example
License
Distributed under the terms of the MIT license, dynkin
is free and open source software
Expand source code
"""
A toolkit for 3D dynamics and kinematics of rigid bodies using the
YPR euler angle convention.
.. include:: ../README.md
"""
__pdoc__ = {
'': False,
'euler': False,
'frame': False,
'rigid_body': False,
}
import numpy as np
def assert_3D_vector(v):
"""
Make sure v is a 3x1 numpy array with dtype==float
"""
out = np.array(v, dtype=float)
if not out.size == 3:
raise AssertionError('Must be a np.ndarray of size 3!')
return out
def assert_6D_vector(v):
"""
Make sure v is a 6x1 numpy array with dtype==float
"""
out = np.array(v, dtype=float)
if not out.size == 6:
raise AssertionError('Must be a np.ndarray of size 6!')
return out
def assert_6D_matrix(m):
"""
Make sure m is a 6x6 numpy array with dtype==float
"""
out = np.array(m, dtype=float)
if not out.shape == (6, 6):
raise ValueError('Must be a np.ndarray of shape (6,6)!')
return out
from .frame import Frame, Transform, transform
from .rigid_body import RigidBody
__all__ = [
'RigidBody',
'Frame',
'Transform',
'transform'
]
Functions
def transform(zeroth_frame, end_frame)
-
Finds a Transform between two frames
Expand source code
def transform(zeroth_frame, end_frame): """ Finds a Transform between two frames """ HTM = np.eye(4) f = end_frame while f: if f == zeroth_frame: return Transform(HTM) HTM = f.HTM @ HTM f = f.parent # If we get to here, the least common base fram is the inertial frame (None) T = Transform(HTM) if zeroth_frame is None: return T T_ = transform(None, zeroth_frame).inv() return Transform(T_.HTM @ T.HTM)
Classes
class Frame (position=None, attitude=None, linear_velocity=None, angular_velocity=None, degrees=False, parent=None)
-
A coordinate reference frame
Expand source code
class Frame: """ A coordinate reference frame """ def __init__(self, position=None, attitude=None, linear_velocity=None, angular_velocity=None, degrees=False, parent=None): self._parent = parent self._HTM = np.eye(4) self.set_position(position or [0, 0, 0]) self.set_attitude(attitude or [0, 0, 0,], degrees=degrees) self.set_linear_velocity(linear_velocity or [0, 0, 0]) self.set_angular_velocity(angular_velocity or [0, 0, 0], degrees=degrees) @property def parent(self): """ Reference to parent Frame """ return self._parent @property def HTM(self): """ Homogenous Transformation Matrix (4x4) defining this Frame """ return self._HTM @property def position(self): """ Position (3x1) of this Frame in relation to parent Frame """ return self._HTM[:3, -1] def set_position(self, position): """ Set position (3x1) of this Frame in relation to parent Frame """ self._HTM[:3, -1] = assert_3D_vector(position) @property def rotation(self): """ Rotation (3x3) of this Frame in relation to parent Frame """ return self.HTM[:3, :3] def set_rotation(self, rotation): """ Set rotation (3x3) of this Frame in relation to parent Frame """ self._HTM[:3, :3] = rotation @property def attitude(self): """ Attitude (euler angles)(3x1) of this Frame in relation to parent Frame """ return rotation_to_euler(self.rotation) def set_attitude(self, attitude, degrees=False): """ Set attitude (euler angles)(3x1) of this Frame in relation to parent Frame """ self._HTM[:3, :3] = euler_to_rotation( assert_3D_vector(np.deg2rad(attitude) if degrees else attitude) ) @property def linear_velocity(self): """ Linear velocity (3x1) of this Frame in relation to parent Frame, decomposed in this Frame """ return self._linear_velocity def set_linear_velocity(self, velocity): """ Set linear velocity (3x1) of this Frame in relation to parent Frame, decomposed in this Frame """ self._linear_velocity = assert_3D_vector(velocity) @property def angular_velocity(self): """ Angular velocity (3x1) of this Frame in relation to parent Frame, decomposed in this Frame """ return self._omega def set_angular_velocity(self, omega, degrees=False): """ Set angular velocity (3x1) of this Frame in relation to parent Frame, decomposed in this Frame """ self._omega = assert_3D_vector(np.deg2rad(omega) if degrees else omega) def get_pose(self): """ Pose (6x1) of this Frame in relation to the inertial Frame """ t = transform(None, self) return np.array([ *t.HTM[:3, -1], *rotation_to_euler(t.HTM[:3, :3]) ]) def get_twist(self): """ Twist (6x1) of this Frame in relation to the inertial Frame """ v = self.linear_velocity w = self.angular_velocity if self.parent is not None: twist_p = self.parent.get_twist() t = transform(self, self.parent) v += t.apply_vector(twist_p[:3] + np.cross(twist_p[3:], self.position)) w += t.apply_vector(twist_p[3:]) return np.array([*v, *w]) def align_child(self, *args, **kwargs): """ Align a child frame in relation to this Frame, accepts same args and kwargs as Frame() """ kwargs.update({ 'parent': self }) return Frame(*args, **kwargs)
Subclasses
- dynkin.rigid_body.RigidBody
Instance variables
var HTM
-
Homogenous Transformation Matrix (4x4) defining this Frame
Expand source code
@property def HTM(self): """ Homogenous Transformation Matrix (4x4) defining this Frame """ return self._HTM
var angular_velocity
-
Angular velocity (3x1) of this Frame in relation to parent Frame, decomposed in this Frame
Expand source code
@property def angular_velocity(self): """ Angular velocity (3x1) of this Frame in relation to parent Frame, decomposed in this Frame """ return self._omega
var attitude
-
Attitude (euler angles)(3x1) of this Frame in relation to parent Frame
Expand source code
@property def attitude(self): """ Attitude (euler angles)(3x1) of this Frame in relation to parent Frame """ return rotation_to_euler(self.rotation)
var linear_velocity
-
Linear velocity (3x1) of this Frame in relation to parent Frame, decomposed in this Frame
Expand source code
@property def linear_velocity(self): """ Linear velocity (3x1) of this Frame in relation to parent Frame, decomposed in this Frame """ return self._linear_velocity
var parent
-
Reference to parent Frame
Expand source code
@property def parent(self): """ Reference to parent Frame """ return self._parent
var position
-
Position (3x1) of this Frame in relation to parent Frame
Expand source code
@property def position(self): """ Position (3x1) of this Frame in relation to parent Frame """ return self._HTM[:3, -1]
var rotation
-
Rotation (3x3) of this Frame in relation to parent Frame
Expand source code
@property def rotation(self): """ Rotation (3x3) of this Frame in relation to parent Frame """ return self.HTM[:3, :3]
Methods
def align_child(self, *args, **kwargs)
-
Align a child frame in relation to this Frame, accepts same args and kwargs as Frame()
Expand source code
def align_child(self, *args, **kwargs): """ Align a child frame in relation to this Frame, accepts same args and kwargs as Frame() """ kwargs.update({ 'parent': self }) return Frame(*args, **kwargs)
def get_pose(self)
-
Pose (6x1) of this Frame in relation to the inertial Frame
Expand source code
def get_pose(self): """ Pose (6x1) of this Frame in relation to the inertial Frame """ t = transform(None, self) return np.array([ *t.HTM[:3, -1], *rotation_to_euler(t.HTM[:3, :3]) ])
def get_twist(self)
-
Twist (6x1) of this Frame in relation to the inertial Frame
Expand source code
def get_twist(self): """ Twist (6x1) of this Frame in relation to the inertial Frame """ v = self.linear_velocity w = self.angular_velocity if self.parent is not None: twist_p = self.parent.get_twist() t = transform(self, self.parent) v += t.apply_vector(twist_p[:3] + np.cross(twist_p[3:], self.position)) w += t.apply_vector(twist_p[3:]) return np.array([*v, *w])
def set_angular_velocity(self, omega, degrees=False)
-
Set angular velocity (3x1) of this Frame in relation to parent Frame, decomposed in this Frame
Expand source code
def set_angular_velocity(self, omega, degrees=False): """ Set angular velocity (3x1) of this Frame in relation to parent Frame, decomposed in this Frame """ self._omega = assert_3D_vector(np.deg2rad(omega) if degrees else omega)
def set_attitude(self, attitude, degrees=False)
-
Set attitude (euler angles)(3x1) of this Frame in relation to parent Frame
Expand source code
def set_attitude(self, attitude, degrees=False): """ Set attitude (euler angles)(3x1) of this Frame in relation to parent Frame """ self._HTM[:3, :3] = euler_to_rotation( assert_3D_vector(np.deg2rad(attitude) if degrees else attitude) )
def set_linear_velocity(self, velocity)
-
Set linear velocity (3x1) of this Frame in relation to parent Frame, decomposed in this Frame
Expand source code
def set_linear_velocity(self, velocity): """ Set linear velocity (3x1) of this Frame in relation to parent Frame, decomposed in this Frame """ self._linear_velocity = assert_3D_vector(velocity)
def set_position(self, position)
-
Set position (3x1) of this Frame in relation to parent Frame
Expand source code
def set_position(self, position): """ Set position (3x1) of this Frame in relation to parent Frame """ self._HTM[:3, -1] = assert_3D_vector(position)
def set_rotation(self, rotation)
-
Set rotation (3x3) of this Frame in relation to parent Frame
Expand source code
def set_rotation(self, rotation): """ Set rotation (3x3) of this Frame in relation to parent Frame """ self._HTM[:3, :3] = rotation
class RigidBody (inertia, cog=None)
-
Idealized rigid body implementation
Expand source code
class RigidBody(Frame): """ Idealized rigid body implementation """ def __init__(self, inertia, cog=None): super().__init__() self._inertia = assert_6D_matrix(inertia) self._cog = self.align_child( position=cog ) @classmethod def from_mass_gyradii(cls, mass, gyradii, cog=None): """ Instantiate from mass and gyradii """ return cls( generalized_inertia_matrix(mass, gyradii), cog=cog ) @property def mass(self): """ Mass of this RigidBody """ return self._inertia[0, 0] @property def CoG(self): """ Reference to Centre of Gravity Frame of this RigidBody """ return self._cog def coriolis_centripetal_matrix(self, inertia, twist): """ Returns Coriolis-Centripetal matrix (6x6) of this RigidBody """ inertia = assert_6D_matrix(inertia) twist = assert_6D_vector(twist) C = np.zeros((6,6)) C[:3, :3] = self.mass*skew(twist[3:]) C[3:, 3:] = -skew(inertia[3:, 3:] @ twist[3:]) return C def generalized_coordinates(self): """ Returns the generalized coordinates of this RigidBody, [x, y, z, fi, theta, psi] """ return self.get_pose() def generalized_velocities(self): """ Returns the generalized velocities of this RigidBody, [dx, dy, dz, dfi, dtheta, dpsi] """ twist = self.get_twist() t = transform(None, self) return np.array([ *t.apply_vector(twist[:3]), *angular_velocity_to_deuler(self.attitude, twist[3:]) ]) def acceleration(self, wrench, additional_inertia=None): """ Returns the resulting acceleration from the given wrench """ f = assert_6D_vector(wrench) additional_inertia = assert_6D_matrix(additional_inertia or np.zeros((6, 6))) twist = self.get_twist() H = motion_transformation_matrix(self.CoG.position) I_cg = self._inertia + additional_inertia C_cg = self.coriolis_centripetal_matrix(I_cg, twist) I_b = H.T @ I_cg @ H C_b = H.T @ C_cg @ H f_cc = C_b @ twist f -= f_cc return np.linalg.solve(I_b, f) def wrench(self, acceleration, additional_inertia=None): """ Returns the required wrench to obtain the given acceleration """ a = assert_6D_vector(acceleration) additional_inertia = assert_6D_matrix(additional_inertia or np.zeros((6, 6))) twist = self.get_twist() H = motion_transformation_matrix(self.CoG.position) I_cg = self._inertia + additional_inertia C_cg = self.coriolis_centripetal_matrix(I_cg, twist) I_b = H.T @ I_cg @ H C_b = H.T @ C_cg @ H f_cc = C_b @ twist f = I_b @ a f += f_cc return f
Ancestors
- dynkin.frame.Frame
Static methods
def from_mass_gyradii(mass, gyradii, cog=None)
-
Instantiate from mass and gyradii
Expand source code
@classmethod def from_mass_gyradii(cls, mass, gyradii, cog=None): """ Instantiate from mass and gyradii """ return cls( generalized_inertia_matrix(mass, gyradii), cog=cog )
Instance variables
var CoG
-
Reference to Centre of Gravity Frame of this RigidBody
Expand source code
@property def CoG(self): """ Reference to Centre of Gravity Frame of this RigidBody """ return self._cog
var mass
-
Mass of this RigidBody
Expand source code
@property def mass(self): """ Mass of this RigidBody """ return self._inertia[0, 0]
Methods
def acceleration(self, wrench, additional_inertia=None)
-
Returns the resulting acceleration from the given wrench
Expand source code
def acceleration(self, wrench, additional_inertia=None): """ Returns the resulting acceleration from the given wrench """ f = assert_6D_vector(wrench) additional_inertia = assert_6D_matrix(additional_inertia or np.zeros((6, 6))) twist = self.get_twist() H = motion_transformation_matrix(self.CoG.position) I_cg = self._inertia + additional_inertia C_cg = self.coriolis_centripetal_matrix(I_cg, twist) I_b = H.T @ I_cg @ H C_b = H.T @ C_cg @ H f_cc = C_b @ twist f -= f_cc return np.linalg.solve(I_b, f)
def coriolis_centripetal_matrix(self, inertia, twist)
-
Returns Coriolis-Centripetal matrix (6x6) of this RigidBody
Expand source code
def coriolis_centripetal_matrix(self, inertia, twist): """ Returns Coriolis-Centripetal matrix (6x6) of this RigidBody """ inertia = assert_6D_matrix(inertia) twist = assert_6D_vector(twist) C = np.zeros((6,6)) C[:3, :3] = self.mass*skew(twist[3:]) C[3:, 3:] = -skew(inertia[3:, 3:] @ twist[3:]) return C
def generalized_coordinates(self)
-
Returns the generalized coordinates of this RigidBody, [x, y, z, fi, theta, psi]
Expand source code
def generalized_coordinates(self): """ Returns the generalized coordinates of this RigidBody, [x, y, z, fi, theta, psi] """ return self.get_pose()
def generalized_velocities(self)
-
Returns the generalized velocities of this RigidBody, [dx, dy, dz, dfi, dtheta, dpsi]
Expand source code
def generalized_velocities(self): """ Returns the generalized velocities of this RigidBody, [dx, dy, dz, dfi, dtheta, dpsi] """ twist = self.get_twist() t = transform(None, self) return np.array([ *t.apply_vector(twist[:3]), *angular_velocity_to_deuler(self.attitude, twist[3:]) ])
def wrench(self, acceleration, additional_inertia=None)
-
Returns the required wrench to obtain the given acceleration
Expand source code
def wrench(self, acceleration, additional_inertia=None): """ Returns the required wrench to obtain the given acceleration """ a = assert_6D_vector(acceleration) additional_inertia = assert_6D_matrix(additional_inertia or np.zeros((6, 6))) twist = self.get_twist() H = motion_transformation_matrix(self.CoG.position) I_cg = self._inertia + additional_inertia C_cg = self.coriolis_centripetal_matrix(I_cg, twist) I_b = H.T @ I_cg @ H C_b = H.T @ C_cg @ H f_cc = C_b @ twist f = I_b @ a f += f_cc return f
class Transform (HTM)
-
A Transform to be applied to positions, vectors, wrenches
Expand source code
class Transform: """ A Transform to be applied to positions, vectors, wrenches """ def __init__(self, HTM): self._HTM = HTM @property def translation(self): """ Translation part (3x1) of Transform """ return self._HTM[:3, -1] @property def HTM(self): """ Homogenous Transformation Matrix (4x4) of Transform """ return self._HTM def apply_vector(self, vector): """ Applies this Transform to vector """ vector = assert_3D_vector(vector) return self._HTM[:3, :3] @ vector def apply_position(self, position): """ Applies this Transform to position """ tmp = np.ones(4) tmp[:3] = assert_3D_vector(position) return (self._HTM @ tmp)[:3] def apply_wrench(self, wrench): """ Applies this Transform to wrench """ wrench = assert_6D_vector(wrench) out = wrench[:] out[:3] = self.apply_vector(wrench[:3]) out[3:] = self.apply_vector(wrench[3:]) + np.cross(self.translation, out[:3]) return out def inv(self): """ Returns inverse of this Transform """ HTM_ = np.linalg.inv(self._HTM) return Transform(HTM_)
Instance variables
var HTM
-
Homogenous Transformation Matrix (4x4) of Transform
Expand source code
@property def HTM(self): """ Homogenous Transformation Matrix (4x4) of Transform """ return self._HTM
var translation
-
Translation part (3x1) of Transform
Expand source code
@property def translation(self): """ Translation part (3x1) of Transform """ return self._HTM[:3, -1]
Methods
def apply_position(self, position)
-
Applies this Transform to position
Expand source code
def apply_position(self, position): """ Applies this Transform to position """ tmp = np.ones(4) tmp[:3] = assert_3D_vector(position) return (self._HTM @ tmp)[:3]
def apply_vector(self, vector)
-
Applies this Transform to vector
Expand source code
def apply_vector(self, vector): """ Applies this Transform to vector """ vector = assert_3D_vector(vector) return self._HTM[:3, :3] @ vector
def apply_wrench(self, wrench)
-
Applies this Transform to wrench
Expand source code
def apply_wrench(self, wrench): """ Applies this Transform to wrench """ wrench = assert_6D_vector(wrench) out = wrench[:] out[:3] = self.apply_vector(wrench[:3]) out[3:] = self.apply_vector(wrench[3:]) + np.cross(self.translation, out[:3]) return out
def inv(self)
-
Returns inverse of this Transform
Expand source code
def inv(self): """ Returns inverse of this Transform """ HTM_ = np.linalg.inv(self._HTM) return Transform(HTM_)