A robot configuration is described by generalized coordinates \( q \in \mathbb{R}^n \), where \( n \) is the number of joint variables and possibly base coordinates. Real mechanisms are rarely completely free: joints, linkages, and contacts impose constraints on admissible motions. Formally, a constraint is a relation among \( q \), \( \dot{q} \), and possibly time \( t \), that every physically realizable motion must satisfy.

At a high level:

• Holonomic constraints restrict configurations and can be written as configuration-level equations \( \phi(q,t)=0 \). They define a lower-dimensional constraint manifold inside the configuration space.
• Non-holonomic constraints restrict velocities and are typically expressed as linear relations in \( \dot{q} \), but they cannot be integrated to equivalent configuration-only equations. They are path-dependent.

This distinction is fundamental for modeling wheeled mobile robots, closed kinematic chains, and robots with rolling or slipping contacts.

flowchart TD
  A["Physical mechanism (joints, contacts, rolling)"] --> B["Mathematical constraint on (q, qdot, t)"]
  B --> C["Can constraint be written as phi(q, t) = 0?"]
  C -->|yes| H["Holonomic constraint"]
  C -->|no| NH["Non-holonomic constraint"]
  H --> H1["Config manifold of lower dimension"]
  NH --> NH1["Velocity distribution; path-dependent reachability"]
        

Let \( \mathcal{Q} \subseteq \mathbb{R}^n \) denote the configuration space of a robot with generalized coordinates \( q = (q_1,\dots,q_n)^\top \). A set of holonomic constraints consists of functions \( \phi_k : \mathcal{Q} \times \mathbb{R} \to \mathbb{R} \), \( k = 1,\dots,m_h \), such that admissible motions satisfy

\[ \phi_k(q(t), t) = 0, \quad k = 1,\dots,m_h. \]

Collecting all constraints into \( \boldsymbol{\phi}(q,t) = [\phi_1(q,t),\dots,\phi_{m_h}(q,t)]^\top \), the admissible configurations at time \( t \) lie on the manifold

\[ \mathcal{Q}_c(t) = \{\, q \in \mathcal{Q} \mid \boldsymbol{\phi}(q,t) = \mathbf{0} \,\}. \]

Assume the Jacobian of constraints with respect to configuration, \( \mathbf{J}_\phi(q,t) = \frac{\partial \boldsymbol{\phi}}{\partial q}(q,t) \), has constant rank \( m_h \) in a neighborhood of interest. Then, by the implicit function theorem, locally we can reparametrize the configuration using \( n - m_h \) independent coordinates, and the number of configuration degrees of freedom is

\[ n_{\text{dof,conf}} = n - m_h. \]

Differentiating the holonomic constraint with respect to time gives its velocity-level counterpart:

\[ \frac{d}{dt}\boldsymbol{\phi}(q(t),t) = \frac{\partial \boldsymbol{\phi}}{\partial q}(q,t)\,\dot{q} + \frac{\partial \boldsymbol{\phi}}{\partial t}(q,t) = \mathbf{0}. \]

This velocity-level relation is not an independent constraint: it is a direct consequence of the configuration-level constraint. Such velocity constraints are called integrable, and the underlying constraints are holonomic.

Example (Point constrained to a circle). A point mass in the plane has configuration \( q = (x,y)^\top \) with unconstrained configuration space \( \mathbb{R}^2 \). Constraining its position to a circle of radius \( R \) centered at the origin yields

\[ \phi(x,y) = x^2 + y^2 - R^2 = 0. \]

The Jacobian is \( \mathbf{J}_\phi(x,y) = [2x \;\; 2y] \), which is rank one for all \( (x,y) \neq (0,0) \). The constrained configuration space is a 1D manifold (the circle), so \( n_{\text{dof,conf}} = 2 - 1 = 1 \). The velocity constraint is

\[ 2x\dot{x} + 2y\dot{y} = 0 \quad \Longleftrightarrow \quad x\dot{x} + y\dot{y} = 0, \]

which is clearly obtained by differentiating \( \phi(x,y)=0 \). Hence this is a holonomic constraint.

Many robotic systems, especially wheeled mobile robots and mechanisms with rolling contacts, are more naturally described by velocity-level constraints that are linear in \( \dot{q} \). A general Pfaffian constraint is of the form

\[ \mathbf{A}(q,t)\,\dot{q} + \mathbf{a}_0(q,t) = \mathbf{0}, \]

where \( \mathbf{A}(q,t) \in \mathbb{R}^{m_v \times n} \) and \( \mathbf{a}_0(q,t) \in \mathbb{R}^{m_v} \). Rows of \( \mathbf{A} \) define linear forms in \( \dot{q} \). At each configuration \( q \), the admissible velocity subspace is

\[ \mathcal{D}(q,t) = \{\, \dot{q} \in \mathbb{R}^n \mid \mathbf{A}(q,t)\,\dot{q} + \mathbf{a}_0(q,t) = \mathbf{0} \,\}. \]

We say that a Pfaffian constraint is holonomic if there exists a set of configuration-level constraint functions \( \boldsymbol{\phi}(q,t) \) such that the Pfaffian constraints are equivalent to \( \frac{d}{dt}\boldsymbol{\phi}(q(t),t)=\mathbf{0} \) along any admissible trajectory. Otherwise, the constraints are non-holonomic.

Thus, non-holonomic constraints cannot be integrated to configuration-only equations. They restrict allowable instantaneous velocities but not the set of reachable configurations. This leads to phenomena like the ability of a car-like robot to reorient itself in a tight parking space despite being unable to move sideways instantaneously.

3.1 Local Integrability Conditions (Single Constraint, 2D Case)

Consider a single Pfaffian constraint in the plane with configuration \( q = (x,y)^\top \):

\[ a_1(x,y)\,\dot{x} + a_2(x,y)\,\dot{y} = 0. \]

Interpreting \( \dot{x} = \frac{dx}{dt} \), \( \dot{y} = \frac{dy}{dt} \), this can be written as a differential 1-form

\[ a_1(x,y)\,dx + a_2(x,y)\,dy = 0. \]

The constraint is holonomic if there exists a scalar function \( \phi(x,y) \) and a nonzero scalar function \( \mu(x,y) \) such that

\[ d\phi(x,y) = \mu(x,y)\,\big( a_1(x,y)\,dx + a_2(x,y)\,dy \big). \]

Writing \( d\phi = \frac{\partial \phi}{\partial x}dx + \frac{\partial \phi}{\partial y}dy \) and comparing coefficients yields

\[ \frac{\partial \phi}{\partial x} = \mu a_1, \quad \frac{\partial \phi}{\partial y} = \mu a_2. \]

A necessary and sufficient condition for the existence of such \( \phi \) and \( \mu \) is

\[ \frac{\partial (\mu a_1)}{\partial y} - \frac{\partial (\mu a_2)}{\partial x} = 0 \]

in a simply connected region. If no such integrating factor \( \mu \) exists, the Pfaffian constraint is non-holonomic.

In higher dimensions, the integrability of a family of Pfaffian constraints is characterized by Frobenius' theorem: roughly speaking, the distribution of admissible velocities \( \mathcal{D}(q) \) must be closed under Lie brackets of vector fields. When that condition fails, the constraint is non-holonomic.

4.1 Holonomic Constraints in Manipulators

Many classical industrial manipulators are modeled as open-chain mechanisms with joint limits. Joint limits such as \( q_i \in [q_i^{\min}, q_i^{\max}] \) are inequality constraints, and they can often be enforced as box constraints in motion planning rather than as equalities in dynamics. However, closed-chain mechanisms (e.g., 4-bar linkages) give rise to holonomic equality constraints of the form

\[ \phi_k(q) = 0, \quad k=1,\dots,m_h, \]

which enforce loop-closure conditions. These constraints reduce the number of independent generalized coordinates that describe the configuration.

4.2 Rolling Disk (Canonical Non-holonomic System)

Consider a disk of radius \( R \) rolling without slipping on a planar surface. A common set of generalized coordinates is \( q = (x, y, \theta, \varphi)^\top \), where:

• \( x, y \): Cartesian position of the contact point on the plane,
• \( \theta \): heading angle of the disk (yaw),
• \( \varphi \): rotation angle of the disk around its axle.

The no-slip, pure-rolling constraints can be expressed as

\[ \dot{x}\cos\theta + \dot{y}\sin\theta - R\dot{\varphi} = 0, \]

\[ \dot{x}\sin\theta - \dot{y}\cos\theta = 0. \]

These are two Pfaffian constraints of the form \( \mathbf{A}(q)\dot{q} = \mathbf{0} \) with \( \mathbf{A}(q) \in \mathbb{R}^{2 \times 4} \). One can show that no configuration-level constraints \( \phi_k(q) = 0 \) exist that generate these velocity constraints as their time derivatives. Intuitively, by exploiting the rolling constraint, the disk can reach a 2D region in \( (x,y) \) despite being unable to move sideways instantaneously. Thus, the constraints are non-holonomic.

4.3 Differential-Drive Robot vs. Omnidirectional Base

A differential-drive robot with configuration \( q = (x,y,\theta)^\top \) and wheel speeds \( \omega_r,\omega_l \) obeys the kinematic model

\[ \dot{x} = v\cos\theta, \quad \dot{y} = v\sin\theta, \quad \dot{\theta} = \omega, \]

where \( v \) and \( \omega \) are functions of the wheel speeds. The non-slip constraint that forbids lateral motion can be written as

\[ -\dot{x}\sin\theta + \dot{y}\cos\theta = 0, \]

which is non-holonomic. In contrast, an omnidirectional base with three or four omni-wheels can realize arbitrary planar velocities \( (\dot{x},\dot{y},\dot{\theta}) \) without such a non-holonomic constraint and is effectively holonomic at the kinematic level.

flowchart TD
  S["Planar robot"] --> D1["Wheeled with lateral no-slip?"]
  D1 -->|yes| NHR["Non-holonomic \n(e.g. diff drive, car)"]
  D1 -->|no| D2["All planar velocities achievable?"]
  D2 -->|yes| HR["Holonomic \n(e.g. omni base)"]
  D2 -->|no| C["Check for other constraints \n(closed chains, guides)"]
        

For a system with \( n \) generalized coordinates and \( m_v \) independent Pfaffian constraints \( \mathbf{A}(q)\dot{q} = \mathbf{0} \) with \( \operatorname{rank}\mathbf{A}(q) = m_v \), the admissible velocity subspace at \( q \) has dimension

\[ \dim \mathcal{D}(q) = n - m_v. \]

For holonomic constraints, this instantaneous dimension coincides with the dimension of the constraint manifold in configuration space. For non-holonomic constraints, however, the configuration space dimension remains \( n \); only the instantaneous velocities are restricted. Over time, the robot can still explore a higher-dimensional set of configurations thanks to the non-commutativity of motion primitives (e.g., forward and rotational motions).

In geometric terms, \( \mathcal{D}(q) \) is a linear subspace of the tangent space \( T_q\mathcal{Q} \). If the distribution \( q \mapsto \mathcal{D}(q) \) is integrable (satisfies Frobenius' condition), the constraints are holonomic; otherwise, they are non-holonomic and often lead to interesting controllability properties.

We now implement the rolling disk constraints in Python and numerically verify the rank of the constraint matrix \( \mathbf{A}(q) \) and the dimension of the admissible velocity subspace. We also integrate a simple kinematic model using these constraints.

import numpy as np

# Rolling disk parameters
R = 0.1  # radius [m]

def A_rolling_disk(q):
    """
    q = [x, y, theta, phi]
    returns A(q) in A(q) qdot = 0
    """
    x, y, theta, phi = q
    A = np.zeros((2, 4))
    # Constraint 1: xdot*cos(theta) + ydot*sin(theta) - R*phidot = 0
    A[0, 0] = np.cos(theta)      # coefficient of xdot
    A[0, 1] = np.sin(theta)      # coefficient of ydot
    A[0, 2] = 0.0               # coefficient of thetadot
    A[0, 3] = -R                # coefficient of phidot

    # Constraint 2: xdot*sin(theta) - ydot*cos(theta) = 0
    A[1, 0] = np.sin(theta)
    A[1, 1] = -np.cos(theta)
    A[1, 2] = 0.0
    A[1, 3] = 0.0
    return A

def admissible_velocity_basis(q):
    """
    Compute a basis for the nullspace of A(q), representing admissible velocities.
    """
    A = A_rolling_disk(q)
    # Use SVD-based nullspace
    U, S, Vt = np.linalg.svd(A)
    rank = np.sum(S > 1e-9)
    null_dim = Vt.shape[0] - rank
    # Last columns of Vt.T span the nullspace
    basis = Vt.T[:, rank:]
    return basis, rank, null_dim

q0 = np.array([0.0, 0.0, 0.0, 0.0])
basis, rankA, dimD = admissible_velocity_basis(q0)
print("rank(A(q0)) =", rankA)
print("dim D(q0)   =", dimD)
print("Basis for admissible velocities at q0:\n", basis)

# Example: constrain motion to "roll forward" with a scalar speed u
def rolling_disk_kinematics(q, u):
    """
    u is the control specifying forward speed along body x-axis.
    """
    x, y, theta, phi = q
    # Forward velocity v = R * phidot
    v = u
    xdot = v * np.cos(theta)
    ydot = v * np.sin(theta)
    thetadot = 0.0
    phidot = v / R
    return np.array([xdot, ydot, thetadot, phidot])

def simulate_rolling_disk(q0, dt=0.01, T=2.0, u=0.2):
    steps = int(T / dt)
    traj = np.zeros((steps + 1, len(q0)))
    traj[0] = q0
    q = q0.copy()
    for k in range(steps):
        qdot = rolling_disk_kinematics(q, u)
        q = q + dt * qdot
        traj[k + 1] = q
    return traj

traj = simulate_rolling_disk(q0)
print("Final configuration after rolling:", traj[-1])
      

In the above code, we:

• Construct the Pfaffian matrix \( \mathbf{A}(q) \) for the rolling disk.
• Compute a basis of the nullspace of \( \mathbf{A}(q) \), which parameterizes admissible velocities.
• Simulate a simple rolling motion using an explicit kinematic model.

In C++, we can use the Eigen library to handle matrices and compute the nullspace of \( \mathbf{A}(q) \). The snippet below focuses on the rolling disk constraints.

#include <iostream>
#include <Eigen/Dense>

using Eigen::MatrixXd;
using Eigen::VectorXd;

MatrixXd A_rolling_disk(const VectorXd& q, double R) {
    // q = [x, y, theta, phi]
    double theta = q(2);
    MatrixXd A(2, 4);
    A.setZero();

    // Constraint 1: xdot*cos(theta) + ydot*sin(theta) - R*phidot = 0
    A(0, 0) = std::cos(theta);
    A(0, 1) = std::sin(theta);
    A(0, 2) = 0.0;
    A(0, 3) = -R;

    // Constraint 2: xdot*sin(theta) - ydot*cos(theta) = 0
    A(1, 0) = std::sin(theta);
    A(1, 1) = -std::cos(theta);
    A(1, 2) = 0.0;
    A(1, 3) = 0.0;

    return A;
}

int main() {
    double R = 0.1;
    VectorXd q(4);
    q << 0.0, 0.0, 0.0, 0.0;

    MatrixXd A = A_rolling_disk(q, R);

    Eigen::JacobiSVD<MatrixXd> svd(A, Eigen::ComputeFullV);
    const auto& S = svd.singularValues();
    MatrixXd V = svd.matrixV();

    int rank = 0;
    double tol = 1e-9;
    for (int i = 0; i < S.size(); ++i) {
        if (S(i) > tol) {
            rank++;
        }
    }
    int null_dim = V.cols() - rank;
    MatrixXd basis = V.block(0, rank, V.rows(), null_dim);

    std::cout << "rank(A)   = " << rank << std::endl;
    std::cout << "null_dim   = " << null_dim << std::endl;
    std::cout << "Nullspace basis:\n" << basis << std::endl;

    return 0;
}
      

This code mirrors the Python implementation and can be integrated into a larger C++ robotics dynamics codebase where \( \mathbf{A}(q) \) is used to project candidate joint velocities onto the admissible subspace.

In Java (without external linear algebra libraries), we can at least implement the constraint evaluation and check whether a given velocity \( \dot{q} \) satisfies the non-holonomic constraint to within a tolerance.

public class RollingDiskConstraints {

    private double R;

    public RollingDiskConstraints(double radius) {
        this.R = radius;
    }

    // q = [x, y, theta, phi]
    // qdot = [xdot, ydot, thetadot, phidot]
    public double[] constraintValues(double[] q, double[] qdot) {
        double theta = q[2];
        double xdot = qdot[0];
        double ydot = qdot[1];
        double phidot = qdot[3];

        double c1 = xdot * Math.cos(theta) + ydot * Math.sin(theta) - R * phidot;
        double c2 = xdot * Math.sin(theta) - ydot * Math.cos(theta);

        return new double[]{c1, c2};
    }

    public boolean isAdmissible(double[] q, double[] qdot, double tol) {
        double[] c = constraintValues(q, qdot);
        return Math.abs(c[0]) < tol && Math.abs(c[1]) < tol;
    }

    public static void main(String[] args) {
        RollingDiskConstraints model = new RollingDiskConstraints(0.1);
        double[] q = new double[]{0.0, 0.0, 0.0, 0.0};
        double[] qdot = new double[]{0.2, 0.0, 0.0, 0.2 / 0.1}; // forward rolling

        boolean ok = model.isAdmissible(q, qdot, 1e-6);
        System.out.println("Is admissible velocity? " + ok);
    }
}
      

This Java code can be extended with matrix utilities or integrated with a Java-based robotics library to analyze non-holonomic constraints in more complex mechanisms.

MATLAB is widely used for prototyping robot kinematics and dynamics. Below is a simple script that sets up the rolling disk kinematics and integrates its motion using ode45. In Simulink, one would typically represent the equations in block form with an integrator and constraint-enforcing blocks.

function rolling_disk_demo()
    R = 0.1;            % radius
    q0 = [0; 0; 0; 0];  % [x; y; theta; phi]
    Tspan = [0 5];
    u = 0.2;            % forward speed along body x-axis

    % Integrate kinematic model
    ode = @(t, q) rolling_disk_kinematics(t, q, u, R);
    [T, Q] = ode45(ode, Tspan, q0);

    figure; plot(Q(:,1), Q(:,2), 'LineWidth', 2);
    xlabel('x [m]'); ylabel('y [m]');
    axis equal; grid on;
    title('Rolling disk trajectory (non-holonomic)');

end

function qdot = rolling_disk_kinematics(~, q, u, R)
    x = q(1); %#ok<NASGU>
    y = q(2); %#ok<NASGU>
    theta = q(3);
    phi = q(4); %#ok<NASGU>

    v = u;
    xdot = v * cos(theta);
    ydot = v * sin(theta);
    thetadot = 0.0;
    phidot = v / R;

    qdot = [xdot; ydot; thetadot; phidot];
end
      

In Simulink, you can:

• Use integrator blocks for \( x,y,\theta,\varphi \).
• Implement the kinematic relationships in a MATLAB Function block or directly with gains and trigonometric blocks.
• Optionally, add assertion blocks to monitor the constraint residuals (e.g. the lateral velocity should remain near zero).

In Mathematica, we can experiment with integrability conditions for a single Pfaffian constraint in 2D. For example, suppose \( a_1(x,y) = y \) and \( a_2(x,y) = -x \), corresponding to the 1-form \( y\,dx - x\,dy = 0 \). The following snippet checks whether there exists an integrating factor \( \mu(x,y) \) such that the integrability condition is satisfied.

(* Define symbolic functions a1, a2 and unknown integrating factor mu[x,y] *)
Clear[x, y, mu];
a1[x_, y_] := y;
a2[x_, y_] := -x;

eq = D[mu[x, y] * a1[x, y], y] - D[mu[x, y] * a2[x, y], x] == 0;

(* Solve PDE for mu; in general this is nontrivial *)
solution = DSolve[eq, mu, {x, y}]

(* If no nontrivial solution exists, the constraint is non-holonomic *)
      

For many choices of \( a_1,a_2 \), including the rolling disk, no globally defined integrating factor exists, confirming the non-holonomic nature of the constraint. Mathematica is particularly useful for analytic explorations of such integrability conditions, complementing the numerical tools of MATLAB and Python.

In this lesson, we formalized the distinction between holonomic and non-holonomic constraints in robotic systems. Holonomic constraints reduce the dimension of the configuration space and can be written as \( \phi(q,t)=0 \), while non-holonomic constraints appear at the velocity level as non-integrable Pfaffian relations \( \mathbf{A}(q,t)\dot{q} + \mathbf{a}_0(q,t) = \mathbf{0} \).

We explored canonical examples such as closed-chain linkages (holonomic) and wheeled mobile robots or rolling disks (non-holonomic), and we introduced the notion of admissible velocity distributions \( \mathcal{D}(q) \subset T_q\mathcal{Q} \). Finally, we implemented these ideas in Python, C++, Java, MATLAB, and Mathematica to highlight how the classification of constraints influences simulation and analysis. In the next lesson, we will incorporate constraints into the Lagrangian formulation using Lagrange multipliers and derive constrained equations of motion.