Introduction

The Yeadon software was developed as a component in our tool chain for performing simulations of bicycle dynamics that include the effects of a human rider. We choose to model the bicycle rider such that they were rigidly attached to the rear frame of the bicycle thus we required estimates of the inertial properties of a human in a seated position. Yeadon's popular human inertia model is a great candidate for this purpose. This notebook describes the process of solving for the joint angles needed to configure a general Yeadon human model to sit on a typical bicycle with hands on the handlebars, arms hanging down, feet at the crank axle, and butt on the seat. Then, we express the inertial properties of the rider in this configuration in terms that allow them to be combined with the inertial properties of the rear frame of the bicycle.

Software Initialization


In [1]:
from IPython.display import Image

import sympy as sym
import sympy.physics.mechanics as me

import numpy as np

from scipy.optimize import fsolve

# This makes sure the mayavi gui isn't blocking in IPython.
%gui wx

from mayavi import mlab

import yeadon


WARNING:traits.has_traits:DEPRECATED: traits.has_traits.wrapped_class, 'the 'implements' class advisor has been deprecated. Use the 'provides' class decorator.

In [2]:
sym.init_printing(print_builtin=False)

In [3]:
%precision 5


Out[3]:
u'%.5f'

Bicycle Geometry

We describe the basic bicycle geometry as follows:

  • $\lambda_{st}$: The acute angle between the seat tube axis and the ground plane
  • $l_{st}$: seat tube length
  • $l_{sp}$: seat post length
  • $l_{cs}$: Projected chain stay length (i.e. projected into the $YZ$ plane of the Yeadon base reference frame)
  • $r_R,r_F$: The rear and front wheel radii.
  • $h_{bb}$: The bottom bracket height.
  • $w_{hb}$: The handlebar width (i.e. measured at the grip point).
  • $w$: wheelbase
  • $L_{hbF}$: The distance from the center of the front wheel to the left handlebar grip.
  • $L_{hbR}$: The distance from the center of the rear wheel to the left handlebar grip.

In [4]:
Image('bicycle-geometry.png')


Out[4]:

Human Geometry and Configuration

The human is described by 95 geometrical measurements using Yeadon's method and we will need a subset of them for the configuration computations:

  • $L_{j3_L}$: distance between knee and hip joint centers
  • $L_{j5_L}$: distance between ankle and hip joint centers
  • $L_{j6_L}$: distance between heel and ankle joint center
  • $L_{s4_L}$: distance between the hip and shoulder levels
  • $L_{s4_w}$: shoulder width
  • $L_{a2_L}$: distance between the elbow and shoulder levels
  • $L_{a4_L}$: distance between the wrist and shoulder levels
  • $L_{a5_L}$: distance between the base of the thumb and wrist
  • $\phi_P$: Somersault angle

In [5]:
Image('http://yeadon.readthedocs.org/en/latest/_images/measurements.png')


Out[5]:

The human's configuration is defined by 21 joint angles, also described in Yeadon's method.


In [6]:
Image('http://yeadon.readthedocs.org/en/latest/_images/configuration.png')


Out[6]:

Problem Statement

Our first objective is to solve for the somersault angle; the leg angles PJ1extension and J1J2flexion; and the arm angles CA1extension, CA1adduction, CA1rotation and A1A2extension that place the rider in the proper configuration on a bicycle. We assume the rider's configuration is symmetric so that the rider's right side configuration is a mirror of the left, for example, PJ1extension is equal to PK1extension.

We first configure the rider by adding geometric points of the bicycle to the main Yeadon reference frame that represent the desired location of the butt, hands, and feet. We set the center of the base of the thumb stadia $L_{a5}$ and $L_{b5}$ to lie on the handlebar grips, the origin of the Yeadon model to be at the top of the seat, and the centers of the heel of the foot stadia $L_{j6}$ and $L_{k6}$ to lie on the bottom bracket axis.

Yeadon reference frame and origin

We first define a reference frame that represents the Yeadon model's base frame. The $\hat{x}$ unit vector is directed from right to left, the $\hat{y}$ unit vector directed anterior to posterior, and the $\hat{z}$ unit vector is directed inferior to superior. This is consistent with the configuration figure above.


In [7]:
yeadon_rf = me.ReferenceFrame('Y')

The origin, $O$, of the Yeadon model is at the base of pelvis.


In [8]:
origin = me.Point('O')

Torso angle

The pelvis, $P$, thorax, $T$, and chest, $C$, segements in Yeadon's model are assumed to be fixed relative to each other. The pelvis is then rotated forward through a somersault angle such that the rider's trunk ($P$,$T$, & $C$) is hunched over with respect to the bicycle. Experimentally, we determine the somersault angle, $\phi_P$, by photographing the seated bicyclist and estimating the angle at which the rider is leaned forward.


In [9]:
somersault = sym.symbols('phi_P')

In [10]:
pelvis = yeadon_rf.orientnew('P', 'Axis', (somersault, yeadon_rf.x))

This rotation and all remaining rotations are body-fixed $X$-$Y$-$Z$ rotations.

Leg angles

To place the heels of the rider, $L_{j6}$ and $L_{k6}$, on the bike pedals, we must solve for the required hip extension and knee flexion angles. The thigh is defined to extend at the hip through the angle PJ1extension and the lower leg to flex at the knee through the angle J1J2flexion. We choose hip adduction to be zero, so the centers of the heels are the same distance from the medial plane as the hip centers are. First, we'll describe the position of the heel from the pelvis. Next, we will describe the position of the bottom bracket from the seat. Then, we will constrain the two positions to be equal.

Location of the heel (using Yeadon parameters)


In [11]:
PJ1extension, J1J2flexion = sym.symbols('phi_J1, phi_J2')

We create two new reference frames for the left thigh, $J1$ and the left lower leg, $J2$.


In [12]:
left_thigh = pelvis.orientnew('J1', 'Axis', (PJ1extension, pelvis.x))
left_lower_leg = left_thigh.orientnew('J2', 'Axis', (J1J2flexion, left_thigh.x))

The heel center is located with respect to the seat via measurements of the subject when standing upright:

  • $L_{j3_L}$: distance between knee and hip joint centers
  • $L_{j5_L}$: distance between ankle and hip joint centers
  • $L_{j6_L}$: distance between heel and ankle joint center

The heel point, $H$, is thuse defined:


In [13]:
thigh_length, Lj5L, Lj6L = sym.symbols('L_j3L, L_j5L, L_j6L')

In [14]:
lower_leg_length = (Lj5L - thigh_length) + Lj6L

In [15]:
heel = origin.locatenew('H', -thigh_length * left_thigh.z - lower_leg_length * left_lower_leg.z)
heel.pos_from(origin)


Out[15]:
$$- L_{j3L}\mathbf{\hat{j1}_z} + (L_{j3L} - L_{j5L} - L_{j6L})\mathbf{\hat{j2}_z}$$

Location of bottom bracket (using bicycle parameters)

The seat point, $S$, is the point at the top of the seat that lies on the seat tube axis, and it is coincident with the origin, $O$, of the Yeadon model.


In [16]:
seat = origin.locatenew('S', 0)

The location from the origin to the bottom bracket can be described by the geometrical measurements:

  • $l_{st}$: seat tube length
  • $l_{sp}$: seat post length
  • $\lambda_{st}$: seat tube angle

In [17]:
seat_tube_length, seat_post_length, seat_tube_angle = sym.symbols('l_st, l_sp, lambda_st')

In [18]:
seat_distance = seat_tube_length + seat_post_length # relative to the bottom bracket
bottom_bracket = me.Point('B')
bottom_bracket.set_pos(origin, -seat_distance * (sym.cos(seat_tube_angle) * yeadon_rf.y + 
                                                 sym.sin(seat_tube_angle) * yeadon_rf.z))

Kinematic constraints for leg angles

We find the hip extension and knee flexion by enforcing that the distance from the heel to the bottom bracket is zero. This gives two nonlinear equations in two unknowns.


In [19]:
leg_zero_vector = bottom_bracket.pos_from(heel)
zero1 = me.dot(leg_zero_vector, yeadon_rf.y)
zero2 = me.dot(leg_zero_vector, yeadon_rf.z)

In [20]:
zero1.trigsimp()


Out[20]:
$$- L_{j3L} \sin{\left (\phi_{J1} + \phi_{P} \right )} + \left(- l_{sp} - l_{st}\right) \cos{\left (\lambda_{st} \right )} + \left(L_{j3L} - L_{j5L} - L_{j6L}\right) \sin{\left (\phi_{J1} + \phi_{J2} + \phi_{P} \right )}$$

In [21]:
zero2.trigsimp()


Out[21]:
$$L_{j3L} \cos{\left (\phi_{J1} + \phi_{P} \right )} + \left(- l_{sp} - l_{st}\right) \sin{\left (\lambda_{st} \right )} + \left(- L_{j3L} + L_{j5L} + L_{j6L}\right) \cos{\left (\phi_{J1} + \phi_{J2} + \phi_{P} \right )}$$

These two equations can be solved for the hip extension angle, $\phi_{J1}$, and the knee flexion angle, $\phi_{J2}$, given all other variables (remember, we measure $\phi_{P}$).

Arm angles

The relationship between the bicycle geometry and the arm's coinfguration angles can be found in a similar manner, although now we require full three dimensional rotations.

Location of the hand (using Yeadon parameters)

The left shoulder, $A1$, is oriented with respect to the trunk through body-fixed $X$-$Y$-$Z$ rotations $\phi_{A1}$, $\theta_{A1}$, and $\psi_{A1}$.


In [22]:
CA1extension, CA1adduction, CA1rotation = sym.symbols('phi_A1, theta_A1, psi_A1')

In [23]:
left_upper_arm = pelvis.orientnew('A1', 'Body',
                                  (CA1extension, CA1adduction, CA1rotation), 'XYZ')

The lower arm is oriented with respect to the upper arm, $A1$, through a simple rotation: A1A2extension or $\phi_{A2}$.


In [24]:
A1A2extension = sym.symbols('phi_A2')

In [25]:
left_lower_arm = left_upper_arm.orientnew('A2', 'Axis', (A1A2extension, left_upper_arm.x))

To define the position of the hand, we require the following measurements, all of which come from the measurements figure above.

  • $L_{s4_L}$: distance between the hip and shoulder levels
  • $L_{s4_w}$: shoulder width
  • $L_{a2_L}$: distance between the elbow and shoulder levels
  • $L_{a4_L}$: distance between the wrist and shoulder levels
  • $L_{a5_L}$: distance between the base of the thumb and wrist

In [26]:
shoulder_width = sym.symbols('L_s4w')
torso_length, upper_arm_length, La4L, La5L = sym.symbols('L_s4L, L_a2L, L_a4L, L_a5L')

In [27]:
shoulder = origin.locatenew('La0', shoulder_width / 2 * pelvis.x + torso_length * pelvis.z)
lower_arm_length = (La4L - upper_arm_length) + La5L
hand = shoulder.locatenew('H', -upper_arm_length * left_upper_arm.z - lower_arm_length * left_lower_arm.z)

Location of handlebar grip (using bicycle parameters)

Now, we locate the handlebar grip points with respect to the seat. But several new bicycle geometric quantities are needed. We choose quantities that are simple to measure in practice, and will use these to solve for quantities that are more natural to use in our model.

  • $l_{cs}$: Projected chain stay length (i.e., projected into the $YZ$ plane of the Yeadon base reference frame).
  • $r_R, r_F$: The rear and front wheel radii.
  • $h_{bb}$: The height of the bottom bracket from the ground.
  • $w_{hb}$: The handlebar width (i.e., measured between the grip points).
  • $w$: wheelbase

We define four unknowns that give the $Y$ and $Z$ measures of the vectors from the front and rear wheel centers to the handlebar grips: $\alpha_{y}$, $\alpha_{z}$, $\beta_{y}$, $\beta_{z}$.


In [28]:
chain_stay_length = sym.symbols('l_cs')
front_wheel_radius, rear_wheel_radius = sym.symbols('r_F, r_R')
bottom_bracket_height = sym.symbols('h_bb')
handle_width = sym.symbols('w_hb')
wheelbase = sym.symbols('w')
alpha_y, alpha_z, beta_y, beta_z = sym.symbols('alpha_y, alpha_z, beta_y, beta_z')

In [29]:
rear_wheel_center = me.Point('cR')
measure_z = rear_wheel_radius - bottom_bracket_height
measure_y = sym.sqrt(chain_stay_length**2 - (measure_z)**2)
rear_wheel_center.set_pos(bottom_bracket,
                          measure_y * yeadon_rf.y +
                          measure_z * yeadon_rf.z)

front_wheel_center = me.Point('cF')
front_wheel_center.set_pos(rear_wheel_center,
                           -rear_wheel_radius * yeadon_rf.z - 
                           wheelbase * yeadon_rf.y +
                           front_wheel_radius * yeadon_rf.z)

grip_defined_from_front = me.Point('hbF')
grip_defined_from_front.set_pos(front_wheel_center,
                                handle_width / 2.0 * yeadon_rf.x +
                                alpha_y * yeadon_rf.y +
                                alpha_z * yeadon_rf.z)

grip_defined_from_rear = me.Point('hbR')
grip_defined_from_rear.set_pos(rear_wheel_center,
                               handle_width / 2.0 * yeadon_rf.x +
                               beta_y * yeadon_rf.y + 
                               beta_z * yeadon_rf.z)

The last two vectors are expressed in terms of unknowns. Now, we'll solve for those unknowns, using measured values for the distances between the wheel centers and the handle grip, $L_{hbF}, L_{hbR}$. We require the seat point we defined in a previous section. We know that the magnitudes of the vectors from the wheel centers to the grip must equate to the measurements.


In [30]:
front_to_grip_distance, rear_to_grip_distance = sym.symbols('L_hbF, L_hbR')

In [31]:
zero3 = grip_defined_from_front.pos_from(front_wheel_center).magnitude() - front_to_grip_distance
zero4 = grip_defined_from_rear.pos_from(rear_wheel_center).magnitude() - rear_to_grip_distance

In [32]:
zero3.trigsimp()


Out[32]:
$$- L_{hbF} + \sqrt{\alpha_{y}^{2} + \alpha_{z}^{2} + 0.25 w_{hb}^{2}}$$

In [33]:
zero4.trigsimp()


Out[33]:
$$- L_{hbR} + \sqrt{\beta_{y}^{2} + \beta_{z}^{2} + 0.25 w_{hb}^{2}}$$

Furthermore, the two handlebar vectors relative to the seat must be equal.


In [34]:
bike_param_zero_vector = (grip_defined_from_front.pos_from(seat) - 
                          grip_defined_from_rear.pos_from(seat))
zero5 = me.dot(bike_param_zero_vector, yeadon_rf.y)
zero6 = me.dot(bike_param_zero_vector, yeadon_rf.z)

In [35]:
zero5


Out[35]:
$$\alpha_{y} - \beta_{y} - w$$

In [36]:
zero6


Out[36]:
$$\alpha_{z} - \beta_{z} + r_{F} - r_{R}$$

These four equations can be solved for $\alpha_{y}$, $\alpha_{z}$, $\beta_{y}$, $\beta_{z}$. Now, we know the position of the handlebar grip point, through grip_defined_from_front or grip_defined_from_rear in terms of the bicycle's geometry.

Kinematic constraints for arm angles

Finally, we generate equations to solve for the arm angles with the additional assumption that the plane defined by axis of the upper arm is always parallel to the ground, $XY$ plane.


In [37]:
arm_zero_vector = hand.pos_from(origin) - grip_defined_from_rear.pos_from(origin)
zero7 = me.dot(arm_zero_vector, yeadon_rf.x)
zero8 = me.dot(arm_zero_vector, yeadon_rf.y)
zero9 = me.dot(arm_zero_vector, yeadon_rf.z)

In [38]:
zero7.trigsimp()


Out[38]:
$$- L_{a2L} \sin{\left (\theta_{A1} \right )} + \frac{L_{s4w}}{2} - 0.5 w_{hb} + \left(\sin{\left (\phi_{A2} \right )} \sin{\left (\psi_{A1} \right )} \cos{\left (\theta_{A1} \right )} + \sin{\left (\theta_{A1} \right )} \cos{\left (\phi_{A2} \right )}\right) \left(L_{a2L} - L_{a4L} - L_{a5L}\right)$$

In [39]:
zero8.trigsimp()


Out[39]:
$$L_{a2L} \sin{\left (\phi_{A1} + \phi_{P} \right )} \cos{\left (\theta_{A1} \right )} - L_{s4L} \sin{\left (\phi_{P} \right )} - \beta_{y} - \sqrt{l_{cs}^{2} - \left(- h_{bb} + r_{R}\right)^{2}} + \left(l_{sp} + l_{st}\right) \cos{\left (\lambda_{st} \right )} + \left(- \left(- \sin{\left (\psi_{A1} \right )} \sin{\left (\theta_{A1} \right )} \sin{\left (\phi_{A1} + \phi_{P} \right )} + \cos{\left (\psi_{A1} \right )} \cos{\left (\phi_{A1} + \phi_{P} \right )}\right) \sin{\left (\phi_{A2} \right )} - \sin{\left (\phi_{A1} + \phi_{P} \right )} \cos{\left (\phi_{A2} \right )} \cos{\left (\theta_{A1} \right )}\right) \left(L_{a2L} - L_{a4L} - L_{a5L}\right)$$

In [40]:
zero9.trigsimp()


Out[40]:
$$- L_{a2L} \cos{\left (\theta_{A1} \right )} \cos{\left (\phi_{A1} + \phi_{P} \right )} + L_{s4L} \cos{\left (\phi_{P} \right )} - \beta_{z} + h_{bb} - r_{R} + \left(l_{sp} + l_{st}\right) \sin{\left (\lambda_{st} \right )} + \left(- \left(\sin{\left (\psi_{A1} \right )} \sin{\left (\theta_{A1} \right )} \cos{\left (\phi_{A1} + \phi_{P} \right )} + \sin{\left (\phi_{A1} + \phi_{P} \right )} \cos{\left (\psi_{A1} \right )}\right) \sin{\left (\phi_{A2} \right )} + \cos{\left (\phi_{A2} \right )} \cos{\left (\theta_{A1} \right )} \cos{\left (\phi_{A1} + \phi_{P} \right )}\right) \left(L_{a2L} - L_{a4L} - L_{a5L}\right)$$

The arm x unit vector is orthogonal to Yeadon z unit vector (i.e. parallel to the ground).


In [41]:
zero10 = me.dot(left_upper_arm.x, yeadon_rf.z)

In [42]:
zero10.trigsimp()


Out[42]:
$$\sin{\left (\psi_{A1} \right )} \sin{\left (\phi_{A1} + \phi_{P} \right )} - \sin{\left (\theta_{A1} \right )} \cos{\left (\psi_{A1} \right )} \cos{\left (\phi_{A1} + \phi_{P} \right )}$$

These four equations allow us to solve for $\phi_{A1}$, $\theta_{A1}$, $\psi_{A}$, and $\phi_{A2}$.

Numerical Solution

We now have 10 equations that we can solve simultaneously for the left leg and left arm angles. This will generally require give good guesses in the numerical solver.


In [43]:
equations = [zero1, zero2, zero3, zero4, zero5, zero6, zero7, zero8, zero9, zero10]

We need to solve the ten equations for $\phi_{J1},\phi_{J2},\phi_{A1},\theta_{A1},\psi_{A1},\phi_{A2},\alpha_{y},\alpha_{z},\beta_{y},\beta_{z}$ given values for all of the measurements. The first two equations for the leg angles are not coupled to the remaining equations for the arms. Also, these are nonlinear equations; we will need to use a solver such as scipy.optimize.fsolve to find the solutions to the equations.

Here we list all the known measurements, as well as all the unknowns for which we are solving. Note that the order is important.


In [44]:
unknown_symbols = (PJ1extension, J1J2flexion, CA1extension, CA1adduction, CA1rotation, A1A2extension,
                   alpha_y, alpha_z, beta_y, beta_z)

bicycle_measurement_symbols = (seat_tube_length, seat_post_length, seat_tube_angle, chain_stay_length, 
                               front_wheel_radius, rear_wheel_radius, bottom_bracket_height, handle_width, 
                               wheelbase, front_to_grip_distance, rear_to_grip_distance)

human_measurement_symbols = (thigh_length, Lj5L, Lj6L, torso_length, shoulder_width, upper_arm_length,
                             La4L, La5L, somersault)

all_symbols = unknown_symbols + bicycle_measurement_symbols + human_measurement_symbols

Now we generate a function that is based on the SymPy expressions and that can be evaluated numerically.


In [45]:
numerical_function = sym.lambdify(all_symbols, equations, modules="numpy")

These are typical values of the measurements, expressed in meters except otherwise noted (from the Browser bicycle: http://moorepants.github.io/dissertation/physicalparameters.html#parameter-tables).


In [46]:
n_lst = 0.53
n_lsp = 0.24
n_lamst = 1.195550538 # radians
n_lcs = 0.46
n_rF = 0.34352982332
n_rR = 0.340958858855
n_hbb = 0.295
n_whb = 0.58
n_w = 1.121
n_LhbF = 0.8930
n_LhbR = 0.9213

In [47]:
args = [n_lst, n_lsp, n_lamst, n_lcs, n_rF, n_rR, n_hbb, n_whb, n_w, n_LhbF, n_LhbR]

These typical measurements in meters from an adult rider for the Yeadon model (Jason's measurements) can be loaded from a YAML file using the Yeadon package.


In [48]:
human = yeadon.Human('JasonYeadonMeas.txt')

We then add the needed measurements to the args variable.


In [49]:
for measurement in ['Lj3L', 'Lj5L', 'Lj6L', 'Ls4L', 'Ls4w', 'La2L', 'La4L', 'La5L']:
    args.append(human.meas[measurement])
#n_somersault = 0.123918377 # radians
n_somersault = 0.175 # try 10 degrees instead of actual measurement because of incompatible geometry
args.append(n_somersault)

Now we must give initial guesses for all of the unknowns. These should be close to the correct solution:


In [50]:
g_PJ1extension = -np.deg2rad(90.0)
g_J1J2flexion = np.deg2rad(75.0)
g_CA1extension = -np.deg2rad(15.0)
g_CA1adduction = np.deg2rad(2.0)
g_CA1rotation = np.deg2rad(2.0)
g_A1A2extension = -np.deg2rad(40.0)
g_alpha_y = n_LhbF * np.cos(np.deg2rad(45.0))
g_alpha_z = n_LhbF * np.sin(np.deg2rad(45.0))
g_beta_y = -n_LhbR * np.cos(np.deg2rad(30.0))
g_beta_z = n_LhbR * np.sin(np.deg2rad(30.0))

guess = [g_PJ1extension, g_J1J2flexion, g_CA1extension, g_CA1adduction, g_CA1rotation, g_A1A2extension,
         g_alpha_y, g_alpha_z, g_beta_y, g_beta_z]

Now we use scipy.optimize.fsolve to hone in on the exact solution.


In [51]:
def my_nonlinear_function(unknowns, knowns):
    """Returns the value of the vector nonlinear function given the measured quantities.
    
    Parameters
    ----------
    unknowns : array_like, shape(10,)
        The numerical values of the unknown quantities, 
        i.e. joint angles and d_ry, d_rz, d_fy, d_fz.
    knowns : array_like, shape(19,)
        The numerical values of the known quantities: human and 
        bicycle geometry.
        
    Returns
    -------
    sol : ndarray, shape(10,)
    
    """
    all_values = np.hstack((unknowns, knowns))
    return np.array(numerical_function(*all_values))

In [52]:
solution = fsolve(my_nonlinear_function, guess, args)

We can see that the solution evaluates the functions to zero.


In [53]:
my_nonlinear_function(solution, args)


Out[53]:
array([ -7.80292e-13,  -1.62204e-12,   3.69660e-12,  -6.42830e-12,
         0.00000e+00,  -5.55112e-17,  -3.93261e-11,   4.25754e-12,
         1.54751e-11,  -1.83789e-10])

In [54]:
for s, v in zip(unknown_symbols, solution):
    print('{} = {:1.2f}'.format(s, v))


phi_J1 = -1.16
phi_J2 = 1.19
phi_A1 = -0.62
theta_A1 = -0.16
psi_A1 = 0.32
phi_A2 = -0.25
alpha_y = 0.54
alpha_z = 0.65
beta_y = -0.58
beta_z = 0.65

The magnitude of the vectors from the front and rear wheel centers to the grip should be equal to the measured distances, $L_{hbf}$ and $L_{hbR}$.


In [55]:
print('{} == {}'.format(np.sqrt((n_whb / 2) ** 2 + solution[-3] ** 2 + solution[-4] ** 2), n_LhbF))


0.893000000004 == 0.893

In [56]:
print('{} == {}'.format(np.sqrt((n_whb / 2) ** 2 + solution[-1] ** 2 + solution[-2] ** 2), n_LhbR))


0.921299999994 == 0.9213

The equations can also be printed as Python code.


In [57]:
for i, eq in enumerate(equations):
    print('zero[{}] = {}'.format(i, sym.printing.lambdarepr.lambdarepr(eq)))
    print('')


zero[0] = L_j3L*(-sin(phi_J1)*cos(phi_P) - sin(phi_P)*cos(phi_J1)) + (-l_sp - l_st)*cos(lambda_st) + (-(-sin(phi_J1)*sin(phi_P) + cos(phi_J1)*cos(phi_P))*sin(phi_J2) + (-sin(phi_J1)*cos(phi_P) - sin(phi_P)*cos(phi_J1))*cos(phi_J2))*(-L_j3L + L_j5L + L_j6L)

zero[1] = L_j3L*(-sin(phi_J1)*sin(phi_P) + cos(phi_J1)*cos(phi_P)) + (-l_sp - l_st)*sin(lambda_st) + ((-sin(phi_J1)*sin(phi_P) + cos(phi_J1)*cos(phi_P))*cos(phi_J2) - (sin(phi_J1)*cos(phi_P) + sin(phi_P)*cos(phi_J1))*sin(phi_J2))*(-L_j3L + L_j5L + L_j6L)

zero[2] = -L_hbF + sqrt(alpha_y**2 + alpha_z**2 + 0.25*w_hb**2)

zero[3] = -L_hbR + sqrt(beta_y**2 + beta_z**2 + 0.25*w_hb**2)

zero[4] = alpha_y - beta_y - w

zero[5] = alpha_z - beta_z + r_F - r_R

zero[6] = -L_a2L*sin(theta_A1) + L_s4w/2 - 0.5*w_hb + (sin(phi_A2)*sin(psi_A1)*cos(theta_A1) + sin(theta_A1)*cos(phi_A2))*(L_a2L - L_a4L - L_a5L)

zero[7] = -L_a2L*(-sin(phi_A1)*cos(phi_P)*cos(theta_A1) - sin(phi_P)*cos(phi_A1)*cos(theta_A1)) - L_s4L*sin(phi_P) - beta_y - sqrt(l_cs**2 - (-h_bb + r_R)**2) - (-l_sp - l_st)*cos(lambda_st) + (-(-(sin(phi_A1)*cos(psi_A1) + sin(psi_A1)*sin(theta_A1)*cos(phi_A1))*sin(phi_P) + (-sin(phi_A1)*sin(psi_A1)*sin(theta_A1) + cos(phi_A1)*cos(psi_A1))*cos(phi_P))*sin(phi_A2) + (-sin(phi_A1)*cos(phi_P)*cos(theta_A1) - sin(phi_P)*cos(phi_A1)*cos(theta_A1))*cos(phi_A2))*(L_a2L - L_a4L - L_a5L)

zero[8] = -L_a2L*(-sin(phi_A1)*sin(phi_P)*cos(theta_A1) + cos(phi_A1)*cos(phi_P)*cos(theta_A1)) + L_s4L*cos(phi_P) - beta_z + h_bb - r_R - (-l_sp - l_st)*sin(lambda_st) + (-((sin(phi_A1)*cos(psi_A1) + sin(psi_A1)*sin(theta_A1)*cos(phi_A1))*cos(phi_P) + (-sin(phi_A1)*sin(psi_A1)*sin(theta_A1) + cos(phi_A1)*cos(psi_A1))*sin(phi_P))*sin(phi_A2) + (-sin(phi_A1)*sin(phi_P)*cos(theta_A1) + cos(phi_A1)*cos(phi_P)*cos(theta_A1))*cos(phi_A2))*(L_a2L - L_a4L - L_a5L)

zero[9] = (sin(phi_A1)*sin(psi_A1) - sin(theta_A1)*cos(phi_A1)*cos(psi_A1))*cos(phi_P) + (sin(phi_A1)*sin(theta_A1)*cos(psi_A1) + sin(psi_A1)*cos(phi_A1))*sin(phi_P)

Configure the Yeadon model

We have completed our first objective! We have oriented the rider to sit on the bicycle. Now, we need to determine the rider's inertia properties? First, we tell the Yeadon human object about the configuration in which we want the rider. The default dictionary can be extracted from the human object with:


In [58]:
cfg_dict = human.CFG.copy()

The left side can be set with the solution:


In [59]:
cfg_dict['PJ1extension'] = solution[0]
cfg_dict['J1J2flexion'] = solution[1]
cfg_dict['CA1extension'] = solution[2]
cfg_dict['CA1adduction'] = solution[3]
cfg_dict['CA1rotation'] = solution[4]
cfg_dict['A1A2extension'] = solution[5]
cfg_dict['somersault'] = n_somersault

And the right side can just be a mirror of the left.


In [60]:
cfg_dict['PK1extension'] = cfg_dict['PJ1extension']
cfg_dict['K1K2flexion'] = cfg_dict['J1J2flexion']
cfg_dict['CB1extension'] = cfg_dict['CA1extension'] 
cfg_dict['CB1abduction'] = -cfg_dict['CA1adduction']
cfg_dict['CB1rotation'] = -cfg_dict['CA1rotation']
cfg_dict['B1B2extension'] = cfg_dict['A1A2extension']

The values of the configuration are:


In [61]:
cfg_dict


Out[61]:
{'A1A2extension': -0.24829391918726648,
 'B1B2extension': -0.24829391918726648,
 'CA1adduction': -0.15756629042399525,
 'CA1extension': -0.61930340001934858,
 'CA1rotation': 0.31840490580450487,
 'CB1abduction': 0.15756629042399525,
 'CB1extension': -0.61930340001934858,
 'CB1rotation': -0.31840490580450487,
 'J1J2flexion': 1.1850079456423874,
 'K1K2flexion': 1.1850079456423874,
 'PJ1adduction': 0.0,
 'PJ1extension': -1.1645088784331983,
 'PK1abduction': 0.0,
 'PK1extension': -1.1645088784331983,
 'PTbending': 0.0,
 'PTsagittalFlexion': 0.0,
 'TCsagittalSpinalFlexion': 0.0,
 'TCspinalTorsion': 0.0,
 'somersault': 0.175,
 'tilt': 0.0,
 'twist': 0.0}

The configuration can the be set with:


In [62]:
human.set_CFG_dict(cfg_dict)

Visualize the configuration

Before visualizing the human in the correct configuration we will add the basic geometry of the bicycle for reference.


In [63]:
def numerical_vector(point, origin, ref_frame):
    """Returns an (3,) length array with the x, y, z coordinates of the
    given point with respect to the given origin in the given reference frame."""
    matrix = point.pos_from(origin).to_matrix(ref_frame)
    f = sym.lambdify(all_symbols, matrix, modules='numpy')
    return np.array(f(*np.hstack((solution, args)))).flatten()

In [64]:
numerical_vec_to_fc = numerical_vector(front_wheel_center, origin, yeadon_rf)
numerical_vec_to_rc = numerical_vector(rear_wheel_center, origin, yeadon_rf)

In [65]:
numerical_vec_to_left_grip = numerical_vector(grip_defined_from_front, origin, yeadon_rf)
numerical_vec_to_bottom_bracket = numerical_vector(bottom_bracket, origin, yeadon_rf)

In [66]:
def add_torus(outer_radius, cross_section_radius, center, orientation):
    """Adds a torus to represent the bicycle wheeels to the scene."""
    torus = mlab.pipeline.parametric_surface()
    torus.function = 'torus'
    torus.parametric_function.ring_radius = outer_radius - cross_section_radius
    torus.parametric_function.cross_section_radius = cross_section_radius
    mlab.pipeline.surface(torus)
    torus.children[0].children[0].actor.actor.orientation = orientation
    torus.children[0].children[0].actor.actor.position = center
    return torus

In [67]:
front_torus = add_torus(n_rF - 0.02, 0.02, numerical_vec_to_fc, np.array([0.0, 90.0, 0.0]))
rear_torus = add_torus(n_rR - 0.02, 0.02, numerical_vec_to_rc, np.array([0.0, 90.0, 0.0]))

Show some important points:


In [68]:
points = np.vstack((numerical_vec_to_left_grip,
                    numerical_vec_to_bottom_bracket,
                    numerical_vec_to_fc,
                    numerical_vec_to_rc))
points


Out[68]:
array([[ 0.29   , -0.40641, -0.01772],
       [ 0.     , -0.28221, -0.71642],
       [ 0.     , -0.94551, -0.66789],
       [ 0.     ,  0.17549, -0.67046]])

In [69]:
null = mlab.points3d(points[:, 0], points[:, 1], points[:, 2],
                     0.125 * np.ones_like(points[:, 0]), scale_factor=1.0)


/home/moorepants/anaconda/envs/yeadon-paper/lib/python2.7/site-packages/traits/has_traits.py:1766: FutureWarning: comparison to `None` will result in an elementwise object comparison in the future.
  setattr( self, name, value )

Draw the Yeadon model:


In [70]:
human.draw(mlabobj=mlab)

Take some orthographic screen shots:


In [71]:
mlab.view(azimuth=0.0, elevation=90.0)
mlab.get_engine().scenes[0].scene.camera.parallel_projection = True
mlab.orientation_axes()
mlab.savefig('side-view.png')
Image('side-view.png') #, width='800px')


Out[71]:

In [72]:
mlab.view(azimuth=-90.0)
mlab.savefig('front-view.png')
Image('front-view.png') #, width='800px')


Out[72]:

Human inertia

The inertial properties of the human in this configuration can be expressed in the Yeadon reference frame with respect to the center of mass of the human:


In [73]:
human.mass


Out[73]:
83.500000000000014

In [74]:
human.center_of_mass


Out[74]:
array([[ 0.     ],
       [-0.14028],
       [ 0.08757]])

In [75]:
human.inertia


Out[75]:
matrix([[ 11.2571 ,   0.     ,   0.     ],
        [  0.     ,  10.98524,  -1.70211],
        [  0.     ,  -1.70211,   2.26922]])

Express inertia in bicycle reference frame

Miejaard et. al 2007 uses the SAE standard for defining the reference frame of a bicycle. The $\hat{x}_b$ unit vector is directed from the back to the front of the bicycle, the $\hat{y}_b$ is directed from the left side of the bicycle to the right side, and the $\hat{z}_b$ is directed from top to bottom. So to rotate a reference frame relative to the Yeadon base reference frame into alignment with the bicycle reference frame, you can first rotate through $\pi$ about $\hat{x}$ and then through $-\frac{\pi}{2}$ about $\hat{z}$. The direction cosine matrix that transforms a vector in blank to a vector in blank is then:


In [76]:
bicycle_rf = yeadon_rf.orientnew('B', 'Space', (sym.pi, -sym.pi / 2, 0), 'XZY')
bicycle_rf.dcm(yeadon_rf)


Out[76]:
$$\left[\begin{matrix}0 & -1 & 0\\-1 & 0 & 0\\0 & 0 & -1\end{matrix}\right]$$

The inertia of the human with respect to the bicycle reference frame and about the human center of mass is then found by:


In [77]:
R = np.array([[0.0, -1.0, 0.0],
              [-1.0, 0.0, 0.0],
              [0.0, 0.0, -1.0]])
human.inertia_transformed(rotmat=R)


Out[77]:
matrix([[ 10.98524,   0.     ,  -1.70211],
        [  0.     ,  11.2571 ,   0.     ],
        [ -1.70211,   0.     ,   2.26922]])

The human center of mass position with respect to the bicycle seat is expressed in the bicycle reference frame is:


In [78]:
R.dot(human.center_of_mass)


Out[78]:
array([[ 0.14028],
       [ 0.     ],
       [-0.08757]])

Combine inertia with bicycle frame

Now if we assume the human is rigidly affixed to the rear bicycle frame we can combine the inertia of the two rigid bodies into a single rigid body by using the parallel axis thereom.

First we form the vector from the bicycle origin (rear wheel contact point) to the Yeadon origin:


In [79]:
bicycle_origin_to_yeadon_origin = origin.pos_from(rear_wheel_center) + rear_wheel_radius * yeadon_rf.z

We can express this vector in the bicycle reference frame:


In [80]:
bicycle_origin_to_yeadon_origin.express(bicycle_rf)


Out[80]:
$$(\sqrt{l_{cs}^{2} - \left(- h_{bb} + r_{R}\right)^{2}} + \left(- l_{sp} - l_{st}\right) \operatorname{cos}\left(\lambda_{st}\right))\mathbf{\hat{b}_x} + (- h_{bb} + \left(- l_{sp} - l_{st}\right) \operatorname{sin}\left(\lambda_{st}\right))\mathbf{\hat{b}_z}$$

And use this expression to compute the location given our bicycle measurements:


In [81]:
compute_yeadon_origin = sym.lambdify(bicycle_measurement_symbols[:7],
                                     bicycle_origin_to_yeadon_origin.to_matrix(bicycle_rf), modules='numpy')

In [82]:
yeadon_origin = np.matrix(compute_yeadon_origin(*args[:7])).reshape(3, 1)
yeadon_origin


Out[82]:
matrix([[ 0.17549],
        [ 0.     ],
        [-1.01142]])

The mass center of the human can then be computed with respect to the bicycle origin and in the bicycle reference frame:


In [83]:
human_mass_center_in_bicycle_rf = yeadon_origin + R.dot(human.center_of_mass)

In [84]:
human_mass_center_in_bicycle_rf


Out[84]:
matrix([[ 0.31577],
        [ 0.     ],
        [-1.09899]])

Now we need some parameters for a the rear frame of the bicycle, inertia, center of mass, and mass. These are taken from the Browser bicycle, see http://moorepants.github.io/dissertation/physicalparameters.html.


In [85]:
n_IBxx = 0.52962890621
n_IBxz = -0.116285607878
n_IByy = 1.3163960125
n_IBzz = 0.756786895402
n_xB = 0.275951285677
n_zB = -0.537842424305
n_mB = 9.86

In [86]:
rear_frame_mass_center = np.matrix([[n_xB],
                                    [0.0],
                                    [n_zB]])
rear_frame_mass_center


Out[86]:
matrix([[ 0.27595],
        [ 0.     ],
        [-0.53784]])

In [87]:
rear_frame_inertia_tensor = np.matrix([[n_IBxx, 0.0, n_IBxz],
                                       [0.0, n_IByy, 0.0],
                                       [n_IBxz, 0.0, n_IBzz]])
rear_frame_inertia_tensor


Out[87]:
matrix([[ 0.52963,  0.     , -0.11629],
        [ 0.     ,  1.3164 ,  0.     ],
        [-0.11629,  0.     ,  0.75679]])

Total mass:


In [88]:
total_mass = n_mB + human.mass
total_mass


Out[88]:
93.360000000000014

The total center of mass of the bicycle rear frame and the rider is:


In [89]:
tmp = human.mass * human_mass_center_in_bicycle_rf
total_center_of_mass = np.matrix([[n_mB * n_xB + tmp[0, 0]],
                                  [n_mB * 0.0 + tmp[1, 0]],
                                  [n_mB * n_zB + tmp[2, 0]]]) / total_mass
total_center_of_mass


Out[89]:
matrix([[ 0.31156],
        [ 0.     ],
        [-1.03972]])

Now use the parallel axis thereom to compute the inertias about the new center of mass and sum them.


In [90]:
tmp = (total_center_of_mass - rear_frame_mass_center)
a = tmp[0, 0]
b = tmp[1, 0]
c = tmp[2, 0]
rear_frame_inertia_about_total_CoM = \
    rear_frame_inertia_tensor + n_mB * np.matrix([[b**2 + c**2, -a * b, -a * c],
                                                  [-a * b, c**2 + a**2, -b * c],
                                                  [-a * c, -b * c, a**2 + b**2]])
rear_frame_inertia_about_total_CoM


Out[90]:
matrix([[ 3.01322,  0.     ,  0.05994],
        [ 0.     ,  3.81249,  0.     ],
        [ 0.05994,  0.     ,  0.76929]])

In [91]:
tmp = (total_center_of_mass - human_mass_center_in_bicycle_rf)
a = tmp[0, 0]
b = tmp[1, 0]
c = tmp[2, 0]
human_inertia_about_total_CoM = \
    human.inertia_transformed(rotmat=R) + human.mass * np.matrix([[b**2 + c**2, -a * b, -a * c],
                                                                  [-a * b, c**2 + a**2, -b * c],
                                                                  [-a * c, -b * c, a**2 + b**2]])
human_inertia_about_total_CoM


Out[91]:
matrix([[ 11.27852,   0.     ,  -1.6813 ],
        [  0.     ,  11.55185,   0.     ],
        [ -1.6813 ,   0.     ,   2.2707 ]])

Finally, we have the combined inertia of the rider and the bicycle rear frame expressed in the bicycle frame about the combined center of mass.


In [92]:
total_inertia = rear_frame_inertia_about_total_CoM + human_inertia_about_total_CoM
total_inertia


Out[92]:
matrix([[ 14.29174,   0.     ,  -1.62136],
        [  0.     ,  15.36434,   0.     ],
        [ -1.62136,   0.     ,   3.03999]])

Software Versions


In [93]:
%install_ext version_information.py


Installed version_information.py. To use it, type:
  %load_ext version_information

In [94]:
%load_ext version_information

In [95]:
%version_information numpy, scipy, sympy, mayavi, yeadon


Out[95]:
SoftwareVersion
Python2.7.9 |Continuum Analytics, Inc.| (default, Dec 15 2014, 10:33:51) [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]
IPython2.3.1
OSposix [linux2]
numpy1.9.1
scipy0.15.1
sympy0.7.6
mayavi4.3.1
yeadon1.2.1
Sun Mar 01 23:24:08 2015 PST