Nonholonomic Constraints¶

Note

You can download this example as a Python script: motion.py or Jupyter Notebook: motion.ipynb.

import sympy as sm
import sympy.physics.mechanics as me
me.init_vprinting(use_latex='mathjax')

class ReferenceFrame(me.ReferenceFrame):

    def __init__(self, *args, **kwargs):

        kwargs.pop('latexs', None)

        lab = args[0].lower()
        tex = r'\hat{{{}}}_{}'

        super(ReferenceFrame, self).__init__(*args,
                                             latexs=(tex.format(lab, 'x'),
                                                     tex.format(lab, 'y'),
                                                     tex.format(lab, 'z')),
                                             **kwargs)
me.ReferenceFrame = ReferenceFrame

Learning Objectives¶

After completing this chapter readers will be able to:

Formulate nonholonomic constraints to constrain motion.
Determine if a nonholonomic constraint is essential and not simply the time derivative of a holonomic constraint.
Formulate a motion constraint for rolling without slip.
Define kinematical differential equations and solve them to put in first order form.
Select different choices of generalized speeds.
Solve for the dependent generalized speeds in terms of the independent generalized speeds.
Calculate the degrees of freedom of a multibody system.

Motion Constraints¶

In Holonomic Constraints, we discussed constraints on the configuration of a system. When defining configuration we are only concerned with the locations of points and how reference frames are oriented. In this chapter, we will consider constraints on the motion of a system. Motion concerns how points and reference frames move in time. Take parallel parking a car as a motivating example, Fig. 27.

_images/motion-parallel.svg — Fig. 27 a) two positions (or configurations) of car 2 relative to cars 1 and 3, b) simplest motion to move car 2 into an empy spot between cars 1 and 3, c) actual motion to move car 2 into the empty spot¶

We know that car 2 can be in either the left or right location in a), i.e. the car’s configuration permits either location. But the motion scenario in b) is not possible. A car cannot move from the left configuration to the right configuration by simply sliding directly to the right (see the note below if you question this). Although, this surely would be nice if we could. A car has wheels and only the front wheels can be steered, so the scenario in c) is a viable motion for the car to end up in the correct final configuration. The car has to move in a specific way to get from one configuration to another. This implies that we have some kind of constraint on the motion but not the configuration. Constraints such as these are called nonholonomic constraints and they take the form:

(100)¶\[\begin{split}\bar{f}_n(\dot{\bar{q}}, \bar{q}, t) = 0 \\ \textrm{ where } \\ \bar{f}_n \in \mathbb{R}^m \\ \bar{q} = \left[ q_1, \ldots, q_n\right]^T \in \mathbb{R}^n\end{split}\]

The \(m\) constraints involve the time derivatives of the generalized coordinates and arise from scalar equations derived from velocities.

Note

We could find a very strong person to push the car sideways, overcoming the very high resisting friction force. It is important to note that any constraint is just a model of a physical phenomena. We know that if we push hard enough and low enough that the car’s lateral motion is not constrained. Also, if the car were on ice, then the nonholomonic constraint may be a poor modeling decision.

Chaplygin Sleigh¶

Take the simple example of the Chaplygin Sleigh, sketched out in Fig. 28. A sleigh can slide along a flat plane, but can only move in the direction it is pointing, much like the wheels of the car above. This system is described by three generalized coordinates \(x,y,\theta\). For the motion to only occur along its body fixed \(\hat{a}_x\) direction, the component of velocity in the body fixed \(\hat{a}_y\) direction must equal zero at all times.

_images/motion-sleigh.svg — Fig. 28 Configuration diagram of a Chaplygin Sleigh. The rectangle \(A\) represents a sleigh moving on a plane. Point \(P\) represents the center of the sleigh.¶

Using SymPy Mechanics we can find the velocity of \(P\) and express it in the \(A\) reference frame:

x, y, theta = me.dynamicsymbols('x, y, theta')

N = me.ReferenceFrame('N')
A = me.ReferenceFrame('A')

A.orient_axis(N, theta, N.z)

O = me.Point('O')
P = me.Point('P')

P.set_pos(O, x*N.x + y*N.y)

O.set_vel(N, 0)

P.vel(N).express(A)

\[\displaystyle (\sin{\left(\theta \right)} \dot{y} + \cos{\left(\theta \right)} \dot{x})\hat{a}_x + (- \sin{\left(\theta \right)} \dot{x} + \cos{\left(\theta \right)} \dot{y})\hat{a}_y\]

The single scalar nonholonomic constraint then takes this form:

(101)¶\[{}^N\bar{v}^P \cdot \hat{a}_y = 0\]

because there can be no velocity component in the \(\hat{a}_y\) direction. With SymPy, this is:

fn = P.vel(N).dot(A.y)
fn

\[\displaystyle - \sin{\left(\theta \right)} \dot{x} + \cos{\left(\theta \right)} \dot{y}\]

How do we know that this is, in fact, a nonholonomic constraint and not simply the time derivative of a holonomic constraint?

Recall one of the four-bar linkage holonomic constraints arising from Eq. (86) and time differentiate it:

t = me.dynamicsymbols._t

q1, q2, q3 = me.dynamicsymbols('q1, q2, q3')
la, lb, lc, ln = sm.symbols('l_a, l_b, l_c, l_n')

fhx = la*sm.cos(q1) + lb*sm.cos(q1 + q2) + lc*sm.cos(q1 + q2 + q3) - ln
sm.trigsimp(fhx.diff(t))

\[\displaystyle - l_{a} \sin{\left(q_{1} \right)} \dot{q}_{1} - l_{b} \left(\dot{q}_{1} + \dot{q}_{2}\right) \sin{\left(q_{1} + q_{2} \right)} - l_{c} \left(\dot{q}_{1} + \dot{q}_{2} + \dot{q}_{3}\right) \sin{\left(q_{1} + q_{2} + q_{3} \right)}\]

This looks like a nonholonomic constraint, i.e. it has time derivatives of the coordinates, but we know that if we integrate this equation with respect to time we can retrieve the original holonomic constraint, so it really isn’t a nonholonomic constraint even though it looks like one.

So if we can integrate \(f_n\) with respect to time and we arrive at a function of only the generalized coordinates and time, then we do not have a nonholonomic constraint, but a holonomic constraint in disguise. Unfortunately, it is not generally possible to integrate \(f_n\) so we must check the integrability of \(f_n\) indirectly.

If \(f_n\) of the sleigh was the time derivative of a holonomic constraint \(f_h\) then it must be able to be expressed in this form:

(102)¶\[f_n = \frac{d f_h}{dt} = \frac{\partial f_h}{\partial x} \frac{dx}{dt} + \frac{\partial f_h}{\partial y} \frac{dy}{dt} + \frac{\partial f_h}{\partial \theta} \frac{d\theta}{dt} + \frac{\partial f_h}{\partial t}\]

and a condition of integrability is that the mixed partial derivatives must commute. By inspection of \(f_n\) we see that we can extract the partial derivatives by collecting the coefficients. SymPy’s coeff() can extract the linear coefficients for us:

dfdx = fn.coeff(x.diff(t))
dfdy = fn.coeff(y.diff(t))
dfdth = fn.coeff(theta.diff(t))

dfdx, dfdy, dfdth

\[\displaystyle \left( - \sin{\left(\theta \right)}, \ \cos{\left(\theta \right)}, \ 0\right)\]

Each pair of mixed partials can be calculated. For example \(\frac{\partial^2 f_h}{\partial y \partial x}\) and \(\frac{\partial^2 f_h}{\partial x \partial y}\):

dfdx.diff(y), dfdy.diff(x)

\[\displaystyle \left( 0, \ 0\right)\]

and the other two pairs:

dfdx.diff(theta), dfdth.diff(x)

\[\displaystyle \left( - \cos{\left(\theta \right)}, \ 0\right)\]

dfdy.diff(theta), dfdth.diff(y)

\[\displaystyle \left( - \sin{\left(\theta \right)}, \ 0\right)\]

We see that to for the last two pairs, the mixed partials do not commute. This proves that \(f_n\) is not integrable and is thus an essential nonholonomic constraint that is not a holonomic constraint in disguise.

Exercise

Check whether the mixed partials of the time derivative of the four-bar linkage constraints commute.

Solution

fnx = fhx.diff(t)
dfdq1 = fnx.diff(q1)
dfdq2 = fnx.diff(q2)
dfdq3 = fnx.diff(q3)

All of the mixed partials are the same:

dfdq1.diff(q2) -  dfdq2.diff(q1)

\[\displaystyle 0\]

dfdq2.diff(q3) - dfdq3.diff(q2)

\[\displaystyle 0\]

dfdq3.diff(q1) - dfdq1.diff(q3)

\[\displaystyle 0\]

All of the mixed partials are the same so this is a holonomic constraint in disguise.

Rolling Without Slip¶

It is quite common to make the modeling assumption that a wheel rolls without slip. A wheel best provides its beneficial properties of rolling and propulsion by ensuring that the friction between the wheel and the surface it rolls on is sufficiently high. This avoids relative motion between a point fixed on the wheel and a point fixed on the surface located at the wheel-surface contact location at any given time. This nature can be modeled by a motion constraint. The key to developing the constraint to ensure there is no relative slip velocity is to identify the correct two points, calculate the velocity of those points, and specify that the relative velocity is zero.

_images/motion-wheel.svg — Fig. 29 A 2D disc \(B\) rolling on a motionless plane \(N\).¶

For example, when a 2D disc \(B\) rolls without slip over a motionless plane \(N\) (Fig. 29), the velocity of a point \(C\) fixed in \(B\) at the contact point with the plane must be zero to ensure no slip when observed from the plane’s reference frame. We can state this mathematically as:

(103)¶\[{}^N\bar{v}^{C} = 0\]

One must be careful about calculating this velocity and recognizing that there are numerous points of possible interest at the same wheel-plane contact location. You may consider these points, for example:

A point \(B_C\) that moves in the plane \(N\) which is always located at the wheel-plane contact location. The coordinate \(q_1\) tracks this point in the figure.
A point \(G_C\) that is fixed in the wheel which follows a cycloid curve as it rolls along.
A point \(G\) that is fixed in the plane which is located at the wheel-plane contact point at any given instance of time.
A point \(C\) that is fixed in the wheel which is located at the wheel-plane contact point at any given instance of time.

A motion constraint that ensures rolling without slip, can only be formed by considering the last two points. The vector constraint equation is:

(104)¶\[{}^N\bar{v}^{C} - {}^N\bar{v}^{G} = 0\]

Point \(G\) is fixed in \(N\) so it has no velocity in \(N\):

(105)¶\[{}^N\bar{v}^{G} = 0\]

Point \(C\) is fixed in \(B\). To determine its velocity, take \(B_o\) to be the wheel center which is also fixed in \(B\). Since both points are fixed in \(B\) we can apply the two point velocity theorem.

(106)¶\[{}^N\bar{v}^{C} = {}^N\bar{v}^{B_o} + {}^N\bar{\omega}^B \times \bar{r}^{C/B_o}\]

We can then use two generalized coordinates to describe the position \(q_1\) (from \(O\) fixed in \(N\)) and rotation \(q_2\) of the wheel. The velocity of the wheel center is then:

(107)¶\[{}^N\bar{v}^{B_o} = \dot{q}_1\hat{n}_x\]

The cross product terms are found with the radius of the wheel with \(r\) and the angular velocity to give the velocity of \(C\):

(108)¶\[\begin{split}{}^N\bar{v}^{C} = & \dot{q}_1\hat{n}_x + \dot{q}_2 \hat{n}_z \times -r\hat{n}_y \\ {}^N\bar{v}^{C} = & \dot{q}_1\hat{n}_x + \dot{q}_2 r \hat{n}_x\end{split}\]

Applying the motion constraint and knowing that \({}^N\bar{v}^{G} = 0\) gives us this scalar constraint equation directly from (108):

(109)¶\[\dot{q}_1 + \dot{q}_2 r = 0\]

This is a scalar constraint equation that ensures rolling without slip and involves the time derivatives of the coordinates. It is integrable and thus actually a holonomic constraint, i.e. \(q_1 + q_2 r = 0\). General rolling without slip in three dimensions will be nonholonomic. Take care to calculate the relative velocities of the two points fixed in each of the bodies in rolling contact that are located at the contact point at that instance of time.

Kinematical Differential Equations¶

In Eq. (100) we show the form of the nonholonomic constraints in terms of \(\dot{\bar{q}}\). Newton’s and Euler’s Second Laws of motion will require calculation of acceleration and angular acceleration respectively. These laws of motion are second order differential equations because it involves second time derivatives of distances and angles. Any second order differential equation can be equivalently represented by two first order differential equations by introducing a new variable for any first derivative terms. We are working towards writing the equations of motion of a multibody system, which will be differential equations that are most useful for simulation when in a first order form. To do this, we now introduce the variables \(\bar{u} = \left[u_1, \ldots, u_n\right]^T\) and define them as linear functions of the time derivatives of the generalized coordinates \(\dot{q}_1, \ldots, \dot{q}_n\). These variables are called generalized speeds. They take the form:

(110)¶\[\bar{u} := \mathbf{Y}_k(\bar{q}, t) \dot{\bar{q}} + \bar{z}_k(\bar{q}, t)\]

\(\bar{u}\) must be chosen such that \(\mathbf{Y}_k\) is invertible. If it is, then we solve for \(\dot{\bar{q}}\) we can write these first order differential equations as such:

(111)¶\[\dot{\bar{q}} = \mathbf{Y}_k^{-1}\left(\bar{u} - \bar{z}_k\right)\]

Eq. (111) are called the kinematical differential equations.

The most common, and always valid, choice of generalized speeds is:

(112)¶\[\bar{u} = \mathbf{I} \dot{\bar{q}}\]

where \(\mathbf{I}\) is the identity matrix. This results in \(u_i = \dot{q}_i\) for \(i=1,\ldots,n\).

Now that we have introduced generalized speeds, the nonholonomic constraints can then be written as:

(113)¶\[\begin{split}\bar{f}_n(\bar{u}, \bar{q}, t) = 0 \\ \textrm{ where } \\ \bar{f}_n \in \mathbb{R}^m \\ \bar{u} = \left[ u_1, \ldots, u_n\right]^T \in \mathbb{R}^n\\ \bar{q} = \left[ q_1, \ldots, q_n\right]^T \in \mathbb{R}^n\end{split}\]

Choosing Generalized Speeds¶

There are many possible choices for generalized speed and you are free to select them as you please, as long as they fit the form of equation (110) and \(\mathbf{Y}_k\) is invertible. Some selections of generalized speeds can reduce the complexity of important velocity expressions and if selected carefully may reduce the complexity of the equations of motion we will derive in a later chapters (see [Mitiguy1996] for examples). To see some examples of selecting generalized speeds, take for example the angular velocity of a reference frame which is oriented with a \(z\textrm{-}x\textrm{-}y\) body fixed orientation:

q1, q2, q3 = me.dynamicsymbols('q1, q2, q3')

A = me.ReferenceFrame('A')
B = me.ReferenceFrame('B')

B.orient_body_fixed(A, (q1, q2, q3), 'ZXY')

A_w_B = B.ang_vel_in(A).simplify()
A_w_B

\[\displaystyle (- \sin{\left(q_{3} \right)} \cos{\left(q_{2} \right)} \dot{q}_{1} + \cos{\left(q_{3} \right)} \dot{q}_{2})\hat{b}_x + (\sin{\left(q_{2} \right)} \dot{q}_{1} + \dot{q}_{3})\hat{b}_y + (\sin{\left(q_{3} \right)} \dot{q}_{2} + \cos{\left(q_{2} \right)} \cos{\left(q_{3} \right)} \dot{q}_{1})\hat{b}_z\]

Choice 1¶

If we choose the simplest definition for the \(u\)’s, i.e. \(u_1=\dot{q}_1\), \(u_2=\dot{q}_2\), and \(u_3=\dot{q}_3\), the angular velocity takes this form:

u1, u2, u3 = me.dynamicsymbols('u1, u2, u3')

t = me.dynamicsymbols._t
qdot = sm.Matrix([q1.diff(t), q2.diff(t), q3.diff(t)])
u = sm.Matrix([u1, u2, u3])

A_w_B = A_w_B.xreplace(dict(zip(qdot, u)))
A_w_B

\[\displaystyle (- u_{1} \sin{\left(q_{3} \right)} \cos{\left(q_{2} \right)} + u_{2} \cos{\left(q_{3} \right)})\hat{b}_x + (u_{1} \sin{\left(q_{2} \right)} + u_{3})\hat{b}_y + (u_{1} \cos{\left(q_{2} \right)} \cos{\left(q_{3} \right)} + u_{2} \sin{\left(q_{3} \right)})\hat{b}_z\]

Yk_plus_zk = qdot
Yk_plus_zk

\[\begin{split}\displaystyle \left[\begin{matrix}\dot{q}_{1}\\\dot{q}_{2}\\\dot{q}_{3}\end{matrix}\right]\end{split}\]

Recall from Solving Linear Systems that the Jacobian is a simple way to extract the coefficients of linear terms into a coefficient matrix for a system of linear equations. In this case, we see that this results in the identity matrix.

Yk = Yk_plus_zk.jacobian(qdot)
Yk

\[\begin{split}\displaystyle \left[\begin{matrix}1 & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right]\end{split}\]

Now find \(\bar{z}_k\) by setting the time derivatives of the generalized coordinates to zero:

qd_zero_repl = dict(zip(qdot, sm.zeros(3, 1)))
qd_zero_repl

\[\displaystyle \left\{ \dot{q}_{1} : 0, \ \dot{q}_{2} : 0, \ \dot{q}_{3} : 0\right\}\]

zk = Yk_plus_zk.xreplace(qd_zero_repl)
zk

\[\begin{split}\displaystyle \left[\begin{matrix}0\\0\\0\end{matrix}\right]\end{split}\]

The linear equation can be solved for the \(\dot{q}\)’s, (Eq. (111)):

sm.Eq(qdot, Yk.LUsolve(u - zk))

\[\begin{split}\displaystyle \left[\begin{matrix}\dot{q}_{1}\\\dot{q}_{2}\\\dot{q}_{3}\end{matrix}\right] = \left[\begin{matrix}u_{1}\\u_{2}\\u_{3}\end{matrix}\right]\end{split}\]

Choice 2¶

Another valid choice is to set the \(u\)’s equal to each measure number of the angular velocity expressed in \(B\):

(114)¶\[\begin{split}u_1 = {}^A\bar{\omega}^B \cdot \hat{b}_x \\ u_2 = {}^A\bar{\omega}^B \cdot \hat{b}_y \\ u_3 = {}^A\bar{\omega}^B \cdot \hat{b}_z\end{split}\]

so that:

(115)¶\[{}^A\bar{\omega}^B = u_1\hat{b}_x + u_2\hat{b}_y + u_3\hat{b}_z\]

A_w_B = B.ang_vel_in(A).simplify()
A_w_B

\[\displaystyle (- \sin{\left(q_{3} \right)} \cos{\left(q_{2} \right)} \dot{q}_{1} + \cos{\left(q_{3} \right)} \dot{q}_{2})\hat{b}_x + (\sin{\left(q_{2} \right)} \dot{q}_{1} + \dot{q}_{3})\hat{b}_y + (\sin{\left(q_{3} \right)} \dot{q}_{2} + \cos{\left(q_{2} \right)} \cos{\left(q_{3} \right)} \dot{q}_{1})\hat{b}_z\]

u1_expr = A_w_B.dot(B.x)
u2_expr = A_w_B.dot(B.y)
u3_expr = A_w_B.dot(B.z)

Yk_plus_zk = sm.Matrix([u1_expr, u2_expr, u3_expr])
Yk_plus_zk

\[\begin{split}\displaystyle \left[\begin{matrix}- \sin{\left(q_{3} \right)} \cos{\left(q_{2} \right)} \dot{q}_{1} + \cos{\left(q_{3} \right)} \dot{q}_{2}\\\sin{\left(q_{2} \right)} \dot{q}_{1} + \dot{q}_{3}\\\sin{\left(q_{3} \right)} \dot{q}_{2} + \cos{\left(q_{2} \right)} \cos{\left(q_{3} \right)} \dot{q}_{1}\end{matrix}\right]\end{split}\]

Yk = Yk_plus_zk.jacobian(qdot)
Yk

\[\begin{split}\displaystyle \left[\begin{matrix}- \sin{\left(q_{3} \right)} \cos{\left(q_{2} \right)} & \cos{\left(q_{3} \right)} & 0\\\sin{\left(q_{2} \right)} & 0 & 1\\\cos{\left(q_{2} \right)} \cos{\left(q_{3} \right)} & \sin{\left(q_{3} \right)} & 0\end{matrix}\right]\end{split}\]

zk = Yk_plus_zk.xreplace(qd_zero_repl)
zk

\[\begin{split}\displaystyle \left[\begin{matrix}0\\0\\0\end{matrix}\right]\end{split}\]

Now we form:

sm.Eq(qdot, sm.trigsimp(Yk.LUsolve(u - zk)))

\[\begin{split}\displaystyle \left[\begin{matrix}\dot{q}_{1}\\\dot{q}_{2}\\\dot{q}_{3}\end{matrix}\right] = \left[\begin{matrix}- \frac{u_{1} \sin{\left(q_{3} \right)} - u_{3} \cos{\left(q_{3} \right)}}{\cos{\left(q_{2} \right)}}\\\frac{u_{1} \cos{\left(2 q_{3} \right)} + u_{1} + u_{3} \sin{\left(2 q_{3} \right)}}{2 \cos{\left(q_{3} \right)}}\\u_{1} \sin{\left(q_{3} \right)} \tan{\left(q_{2} \right)} + u_{2} - u_{3} \cos{\left(q_{3} \right)} \tan{\left(q_{2} \right)}\end{matrix}\right]\end{split}\]

Note

Notice how the kinematical differential equations are not valid when \(q_2\) or \(q_3\) are even multiples of \(\pi/2\). If your system must orient through these values, you’ll need to select a different body fixed rotation or an orientation method that isn’t suseptible to these issues.

Choice 3¶

Another valid choice is to set the \(u\)’s equal to each measure number of the angular velocity expressed in \(A\):

(116)¶\[\begin{split}u_1 = {}^A\bar{\omega}^B \cdot \hat{a}_x \\ u_2 = {}^A\bar{\omega}^B \cdot \hat{a}_y \\ u_3 = {}^A\bar{\omega}^B \cdot \hat{a}_z\end{split}\]

so that:

(117)¶\[{}^A\bar{\omega}^B = u_1\hat{a}_x + u_2\hat{a}_y + u_3\hat{a}_z\]

A_w_B = B.ang_vel_in(A).express(A).simplify()
A_w_B

\[\displaystyle (- \sin{\left(q_{1} \right)} \cos{\left(q_{2} \right)} \dot{q}_{3} + \cos{\left(q_{1} \right)} \dot{q}_{2})\hat{a}_x + (\sin{\left(q_{1} \right)} \dot{q}_{2} + \cos{\left(q_{1} \right)} \cos{\left(q_{2} \right)} \dot{q}_{3})\hat{a}_y + (\sin{\left(q_{2} \right)} \dot{q}_{3} + \dot{q}_{1})\hat{a}_z\]

u1_expr = A_w_B.dot(A.x)
u2_expr = A_w_B.dot(A.y)
u3_expr = A_w_B.dot(A.z)

Yk_plus_zk = sm.Matrix([u1_expr, u2_expr, u3_expr])
Yk_plus_zk

\[\begin{split}\displaystyle \left[\begin{matrix}- \sin{\left(q_{1} \right)} \cos{\left(q_{2} \right)} \dot{q}_{3} + \cos{\left(q_{1} \right)} \dot{q}_{2}\\\sin{\left(q_{1} \right)} \dot{q}_{2} + \cos{\left(q_{1} \right)} \cos{\left(q_{2} \right)} \dot{q}_{3}\\\sin{\left(q_{2} \right)} \dot{q}_{3} + \dot{q}_{1}\end{matrix}\right]\end{split}\]

Yk = Yk_plus_zk.jacobian(qdot)
Yk

\[\begin{split}\displaystyle \left[\begin{matrix}0 & \cos{\left(q_{1} \right)} & - \sin{\left(q_{1} \right)} \cos{\left(q_{2} \right)}\\0 & \sin{\left(q_{1} \right)} & \cos{\left(q_{1} \right)} \cos{\left(q_{2} \right)}\\1 & 0 & \sin{\left(q_{2} \right)}\end{matrix}\right]\end{split}\]

zk = Yk_plus_zk.xreplace(qd_zero_repl)
zk

\[\begin{split}\displaystyle \left[\begin{matrix}0\\0\\0\end{matrix}\right]\end{split}\]

sm.Eq(qdot, sm.trigsimp(Yk.LUsolve(u - zk)))

\[\begin{split}\displaystyle \left[\begin{matrix}\dot{q}_{1}\\\dot{q}_{2}\\\dot{q}_{3}\end{matrix}\right] = \left[\begin{matrix}\left(u_{1} \sin{\left(q_{1} \right)} - u_{2} \cos{\left(q_{1} \right)}\right) \tan{\left(q_{2} \right)} + u_{3}\\u_{1} \cos{\left(q_{1} \right)} + u_{2} \sin{\left(q_{1} \right)}\\- \frac{u_{1} \sin{\left(q_{1} \right)} - u_{2} \cos{\left(q_{1} \right)}}{\cos{\left(q_{2} \right)}}\end{matrix}\right]\end{split}\]

Snakeboard¶

A snakeboard is a variation on a skateboard that can be propelled via nonholonomic locomotion [Ostrowski1994]. Similar to the Chaplygin Sleigh, the wheels can generally only travel in the direction they are pointed. This classic video from 1993 shows how to propel the board:

Fig. 30 shows what a real Snakeboard looks like and Fig. 31 shows a configuration diagram.

_images/motion-snakeboard.svg — Fig. 31 Configuration diagram of a planar Snakeboard model.¶

Start by defining the time varying variables and constants:

q1, q2, q3, q4, q5 = me.dynamicsymbols('q1, q2, q3, q4, q5')
l = sm.symbols('l')

The reference frames are all simple rotations about the axis normal to the plane:

N = me.ReferenceFrame('N')
A = me.ReferenceFrame('A')
B = me.ReferenceFrame('B')
C = me.ReferenceFrame('C')

A.orient_axis(N, q3, N.z)
B.orient_axis(A, q4, A.z)
C.orient_axis(A, q5, A.z)

The angular velocities of each reference frame are then:

A.ang_vel_in(N)

\[\displaystyle \dot{q}_{3}\hat{n}_z\]

B.ang_vel_in(N)

\[\displaystyle \dot{q}_{4}\hat{a}_z + \dot{q}_{3}\hat{n}_z\]

C.ang_vel_in(N)

\[\displaystyle \dot{q}_{5}\hat{a}_z + \dot{q}_{3}\hat{n}_z\]

Establish the position vectors among the points:

O = me.Point('O')
Ao = me.Point('A_o')
Bo = me.Point('B_o')
Co = me.Point('C_o')

Ao.set_pos(O, q1*N.x + q2*N.y)
Bo.set_pos(Ao, l/2*A.x)
Co.set_pos(Ao, -l/2*A.x)

The velocity of \(A_o\) in \(N\) is a simple time derivative:

O.set_vel(N, 0)
Ao.vel(N)

\[\displaystyle \dot{q}_{1}\hat{n}_x + \dot{q}_{2}\hat{n}_y\]

The two point theorem is handy for computing the other two velocities:

Bo.v2pt_theory(Ao, N, A)

\[\displaystyle \dot{q}_{1}\hat{n}_x + \dot{q}_{2}\hat{n}_y + \frac{l \dot{q}_{3}}{2}\hat{a}_y\]

Co.v2pt_theory(Ao, N, A)

\[\displaystyle \dot{q}_{1}\hat{n}_x + \dot{q}_{2}\hat{n}_y - \frac{l \dot{q}_{3}}{2}\hat{a}_y\]

The unit vectors of \(B\) and \(C\) are aligned with the wheels of the Snakeboard. This lets us impose that there is no velocity in the direction normal to the wheel’s rolling direction by taking dot products with the respectively reference frames’ \(y\) direction unit vector to form the two nonholonomic constraints:

(118)¶\[\begin{split}{}^A\bar{v}^{Bo} \cdot \hat{b}_y = 0 \\ {}^A\bar{v}^{Co} \cdot \hat{c}_y = 0\end{split}\]

fn = sm.Matrix([Bo.vel(N).dot(B.y),
                Co.vel(N).dot(C.y)])
fn = sm.trigsimp(fn)
fn

\[\begin{split}\displaystyle \left[\begin{matrix}\frac{l \cos{\left(q_{4} \right)} \dot{q}_{3}}{2} - \sin{\left(q_{3} + q_{4} \right)} \dot{q}_{1} + \cos{\left(q_{3} + q_{4} \right)} \dot{q}_{2}\\- \frac{l \cos{\left(q_{5} \right)} \dot{q}_{3}}{2} - \sin{\left(q_{3} + q_{5} \right)} \dot{q}_{1} + \cos{\left(q_{3} + q_{5} \right)} \dot{q}_{2}\end{matrix}\right]\end{split}\]

Now we introduce some generalized speeds. By inspection of the above constraint equations, we can see that defining a generalized speed equal to \(\frac{l\dot{q}_3}{2}\) can simplify the equations a bit. So define these generalized speeds:

(119)¶\[\begin{split}u_i = \dot{q}_i \textrm{ for } i=1,2,4,5 \\ u_3 = \frac{l\dot{q}_3}{2}\end{split}\]

Now replace all of the time derivatives of the generalized coordinates with the generalized speeds. We use subs() here because the replacement isn’t an exact replacement (in the sense of xreplace()).

u1, u2, u3, u4, u5 = me.dynamicsymbols('u1, u2, u3, u4, u5')

u_repl = {
    q1.diff(): u1,
    q2.diff(): u2,
    l*q3.diff()/2: u3,
    q4.diff(): u4,
    q5.diff(): u5
}

fn = fn.subs(u_repl)
fn

\[\begin{split}\displaystyle \left[\begin{matrix}- u_{1} \sin{\left(q_{3} + q_{4} \right)} + u_{2} \cos{\left(q_{3} + q_{4} \right)} + u_{3} \cos{\left(q_{4} \right)}\\- u_{1} \sin{\left(q_{3} + q_{5} \right)} + u_{2} \cos{\left(q_{3} + q_{5} \right)} - u_{3} \cos{\left(q_{5} \right)}\end{matrix}\right]\end{split}\]

These nonholonomic constraints take this form:

(120)¶\[\bar{f}_n(u_1, u_2, u_3, q_3, q_4, q_5) = 0 \textrm{ where } \bar{f}_n \in \mathbb{R}^2\]

We now have two equations with three unknown generalized speeds. Note that all of the generalized coordinates are not present in these constraints which is common. We can solve for two of the generalized speeds in terms of the third. So we select two as dependent generalized speeds and one as an independent generalized speed. Because nonholonomic constraints are derived from measure numbers of velocity vectors, the nonholonomic constraints are always linear in the generalized speeds. If we introduce \(\bar{u}_s\) as a vector of independent generalized speeds and \(\bar{u}_r\) as a vector of dependent generalized speeds, the nonholonomic constraints can be written as:

(121)¶\[\bar{f}_n(\bar{u}_s, \bar{u}_r, \bar{q}, t) = \mathbf{A}_r(\bar{q}, t) \bar{u}_r + \mathbf{A}_s(\bar{q}, t) \bar{u}_s + \bar{b}_{rs}(\bar{q}, t) = 0\]

or

(122)¶\[\begin{split}\bar{u}_r = \mathbf{A}_r^{-1}\left(-\mathbf{A}_s \bar{u}_s - \bar{b}_{rs}\right) \\ \bar{u}_r = \mathbf{A}_n \bar{u}_s + \bar{b}_n\end{split}\]

For the Snakeboard let’s choose \(\bar{u}_s = [u_3, u_4, u_5]^T\) as the independent generalized speeds and \(\bar{u}_r = [u_1, u_2]^T\) as the dependent generalized speeds.

us = sm.Matrix([u3, u4, u5])
ur = sm.Matrix([u1, u2])

\(\mathbf{A}_r\) are the linear coefficients of \(\bar{u}_r\) so:

Ar = fn.jacobian(ur)
Ar

\[\begin{split}\displaystyle \left[\begin{matrix}- \sin{\left(q_{3} + q_{4} \right)} & \cos{\left(q_{3} + q_{4} \right)}\\- \sin{\left(q_{3} + q_{5} \right)} & \cos{\left(q_{3} + q_{5} \right)}\end{matrix}\right]\end{split}\]

\(\mathbf{A}_s\) are the linear coefficients of \(\bar{u}_s\) so:

As = fn.jacobian(us)
As

\[\begin{split}\displaystyle \left[\begin{matrix}\cos{\left(q_{4} \right)} & 0 & 0\\- \cos{\left(q_{5} \right)} & 0 & 0\end{matrix}\right]\end{split}\]

\(\bar{b}_{rs}\) remains when \(\bar{u}=0\):

brs = fn.xreplace(dict(zip([u1, u2, u3, u4, u5], [0, 0, 0, 0, 0])))
brs

\[\begin{split}\displaystyle \left[\begin{matrix}0\\0\end{matrix}\right]\end{split}\]

\(\mathbf{A}_n\) and \(\bar{b}_n\) are formed by solving the linear system:

An = Ar.LUsolve(-As)
An = sm.simplify(An)
An

\[\begin{split}\displaystyle \left[\begin{matrix}\frac{\frac{\cos{\left(q_{3} - q_{4} + q_{5} \right)}}{2} + \frac{\cos{\left(q_{3} + q_{4} - q_{5} \right)}}{2} + \cos{\left(q_{3} + q_{4} + q_{5} \right)}}{\sin{\left(q_{4} - q_{5} \right)}} & 0 & 0\\\frac{\frac{\sin{\left(q_{3} - q_{4} + q_{5} \right)}}{2} + \frac{\sin{\left(q_{3} + q_{4} - q_{5} \right)}}{2} + \sin{\left(q_{3} + q_{4} + q_{5} \right)}}{\sin{\left(q_{4} - q_{5} \right)}} & 0 & 0\end{matrix}\right]\end{split}\]

bn = Ar.LUsolve(-brs)
bn

\[\begin{split}\displaystyle \left[\begin{matrix}0\\0\end{matrix}\right]\end{split}\]

We now have the \(m=2\) dependent generalized speeds \(\bar{u}_r=\left[u_1,u_2\right]^T\) written as functions of the n=1 independent generalized speeds \(\bar{u}_s=\left[u_3\right]\):

sm.Eq(ur, An*us + bn)

\[\begin{split}\displaystyle \left[\begin{matrix}u_{1}\\u_{2}\end{matrix}\right] = \left[\begin{matrix}\frac{\left(\frac{\cos{\left(q_{3} - q_{4} + q_{5} \right)}}{2} + \frac{\cos{\left(q_{3} + q_{4} - q_{5} \right)}}{2} + \cos{\left(q_{3} + q_{4} + q_{5} \right)}\right) u_{3}}{\sin{\left(q_{4} - q_{5} \right)}}\\\frac{\left(\frac{\sin{\left(q_{3} - q_{4} + q_{5} \right)}}{2} + \frac{\sin{\left(q_{3} + q_{4} - q_{5} \right)}}{2} + \sin{\left(q_{3} + q_{4} + q_{5} \right)}\right) u_{3}}{\sin{\left(q_{4} - q_{5} \right)}}\end{matrix}\right]\end{split}\]

Degrees of Freedom¶

For simple nonholonomic systems observed in a reference frame \(A\), such as the Chaplygin Sleigh or the Snakeboard, the degrees of freedom in \(A\) are equal to the number of independent generalized speeds. The number of degrees of freedom \(p\) is defined as:

(123)¶\[p := n - m\]

where \(n\) is the number of generalized coordinates and \(m\) are the number of nonholonomic constraints (and thus dependent generalized speeds). If there are no nonholonomic constraints, the system is a holonomic system in \(A\) and \(p=n\) making the number of degrees of freedom equal to the number of generalized coordinates.

Exercise

What are the number of degrees of freedom for the Chaplygin Sleigh, Snakeboard, and Four-bar linkage?

Solution

The Chapylgin Sleigh has \(n=3\) generalized coordinates \(x,y,\theta\) and \(m=1\) nonholonomic constraints. The degrees of freedom are then \(p = 3 - 1 = 2\).

The Snakeboard has \(n=5\) generalized coordinates and \(m=2\) nonholonomic constraints. The degrees of freedom are then \(p = 5 - 2 = 3\).

We described the four-bar linkage with \(N=3\) coordinates and there were \(M=2\) holonomic constraints leaving us with \(n=N-M=3-2=1\) generalized coordinates. There are no nonholonomic constraints so \(m=0\). This means that there \(p=n-m=1-0=1\) degrees of freedom.

It is not always easy to visualize the degrees of freedom of a nonholonomic system when thinking of its motion, but for holonomic systems thought experiments where you vary one or two generalized coordinates at a time can help you visualize the motion.

If you have a holonomic system (no nonholonomic constraints) the degrees of freedom are equal to the number of generalized coordinates. Nonholonomic systems (those with non-integrable motion constraints) have fewer degrees of freedom than the number of generalized coordinates.