Consider a mechanical system with generalized coordinates \( \mathbf{q}\in\mathbb{R}^n \) and Lagrangian \( L(\mathbf{q},\dot{\mathbf{q}},t) = T(\mathbf{q},\dot{\mathbf{q}}) - V(\mathbf{q}) \). Without constraints, the Lagrange–Euler equations read

\[ \frac{\mathrm{d}}{\mathrm{d}t} \left( \frac{\partial L}{\partial \dot{\mathbf{q}}} \right) - \frac{\partial L}{\partial \mathbf{q}} = \mathbf{Q}, \]

where \( \mathbf{Q} \) collects generalized nonconservative forces (including joint torques \tau in manipulators).

Suppose now that the motion of the system is restricted by holonomic constraints \( \boldsymbol{\phi}(\mathbf{q},t)=\mathbf{0} \), with \( \boldsymbol{\phi}:\mathbb{R}^n\times\mathbb{R}\to\mathbb{R}^m \) and \( m \leq n \). These constraints define a lower-dimensional configuration manifold on which the motion must lie. The constraint forces required to keep the system on this manifold are a priori unknown and are not expressible as gradients of a scalar potential. Lagrange multipliers provide a systematic way to:

• Represent these constraint forces as generalized forces \( \mathbf{Q}^{\mathrm{c}} = \mathbf{J}_\phi^\top(\mathbf{q},t)\boldsymbol{\lambda} \),
• Ensure that the configuration constraints and their time derivatives are satisfied,
• Augment the robot dynamics to solve simultaneously for accelerations \( \ddot{\mathbf{q}} \) and multipliers \( \boldsymbol{\lambda} \).

flowchart TD
  Q["Choose generalized coordinates q"] --> PHI["Specify constraints phi(q,t)=0"]
  PHI --> L["Form Lagrangian L = T - V"]
  L --> EU["Write unconstrained Lagrange-Euler equations"]
  EU --> AUG["Add constraint forces J_phi^T * lambda"]
  AUG --> SYS["Obtain augmented system for qddot and lambda"]
  SYS --> INT["Integrate state (q, qdot) over time"]
        

We recall d'Alembert’s principle for an unconstrained system: for any virtual displacement \( \delta\mathbf{q} \),

\[ \sum_{i=1}^n \left( Q_i - \frac{\mathrm{d}}{\mathrm{d}t} \frac{\partial T}{\partial \dot{q}_i} + \frac{\partial T}{\partial q_i} + \frac{\partial V}{\partial q_i} \right) \delta q_i = 0. \]

For an ideal holonomically constrained system with \( \boldsymbol{\phi}(\mathbf{q},t)=\mathbf{0} \), the virtual displacements must satisfy the virtual constraint

\[ \mathbf{J}_\phi(\mathbf{q},t)\,\delta\mathbf{q} = \mathbf{0},\quad \mathbf{J}_\phi(\mathbf{q},t) := \frac{\partial \boldsymbol{\phi}}{\partial \mathbf{q}}(\mathbf{q},t) \in\mathbb{R}^{m\times n}, \]

where \( \mathbf{J}_\phi \) is the constraint Jacobian. Let \( \mathbf{Q}^{\mathrm{c}} \) denote the generalized constraint forces. d'Alembert’s principle with constraints becomes

\[ \sum_{i=1}^n \left( Q_i + Q_i^{\mathrm{c}} - \frac{\mathrm{d}}{\mathrm{d}t} \frac{\partial T}{\partial \dot{q}_i} + \frac{\partial T}{\partial q_i} + \frac{\partial V}{\partial q_i} \right) \delta q_i = 0 \quad \text{for all } \delta\mathbf{q} \text{ with } \mathbf{J}_\phi\,\delta\mathbf{q} = \mathbf{0}. \]

Because the virtual displacements are restricted by \( \mathbf{J}_\phi\,\delta\mathbf{q} = \mathbf{0} \), the only way the virtual work of the constraint forces can vanish for all admissible \( \delta\mathbf{q} \) is if \( \mathbf{Q}^{\mathrm{c}} \) lies in the span of the rows of \( \mathbf{J}_\phi \):

\[ \mathbf{Q}^{\mathrm{c}}(\mathbf{q},\dot{\mathbf{q}},t) = \mathbf{J}_\phi(\mathbf{q},t)^\top \boldsymbol{\lambda}(\mathbf{q},\dot{\mathbf{q}},t) \]

for some multiplier vector \( \boldsymbol{\lambda}\in\mathbb{R}^m \). Substituting into d'Alembert’s principle and using the usual equivalence between vanishing virtual work and Euler–Lagrange equations, we obtain

\[ \frac{\mathrm{d}}{\mathrm{d}t} \left( \frac{\partial L}{\partial \dot{\mathbf{q}}} \right) - \frac{\partial L}{\partial \mathbf{q}} = \mathbf{Q} + \mathbf{J}_\phi(\mathbf{q},t)^\top\boldsymbol{\lambda}, \quad \boldsymbol{\phi}(\mathbf{q},t)=\mathbf{0}. \]

The multipliers \( \boldsymbol{\lambda} \) encode the constraint reaction forces, while the constraints themselves are enforced explicitly through \( \boldsymbol{\phi}(\mathbf{q},t)=\mathbf{0} \).

For an n-DOF manipulator with generalized coordinates \( \mathbf{q}\in\mathbb{R}^n \), we have already seen the unconstrained dynamics in standard form

\[ \mathbf{M}(\mathbf{q})\,\ddot{\mathbf{q}} + \mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\,\dot{\mathbf{q}} + \mathbf{g}(\mathbf{q}) = \mathbf{S}^\top \tau, \]

where \( \mathbf{M}(\mathbf{q}) \) is the positive-definite inertia matrix, \( \mathbf{C}(\mathbf{q},\dot{\mathbf{q}}) \) is the Coriolis/centrifugal matrix, \( \mathbf{g}(\mathbf{q}) \) is gravity, and \( \mathbf{S} \) selects actuated joints. Now suppose the manipulator is subject to m holonomic constraints

\[ \boldsymbol{\phi}(\mathbf{q},t)=\mathbf{0},\quad \mathbf{J}_\phi(\mathbf{q},t) = \frac{\partial \boldsymbol{\phi}}{\partial \mathbf{q}} \in\mathbb{R}^{m\times n}. \]

The constrained dynamics become

\[ \mathbf{M}(\mathbf{q})\,\ddot{\mathbf{q}} + \mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\,\dot{\mathbf{q}} + \mathbf{g}(\mathbf{q}) = \mathbf{S}^\top \tau + \mathbf{J}_\phi(\mathbf{q},t)^\top\boldsymbol{\lambda}, \quad \boldsymbol{\phi}(\mathbf{q},t)=\mathbf{0}. \]

Differentiating the constraints yields velocity and acceleration-level relations. First derivative:

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

Differentiating once more:

\[ \mathbf{J}_\phi(\mathbf{q},t)\,\ddot{\mathbf{q}} + \dot{\mathbf{J}}_\phi(\mathbf{q},\dot{\mathbf{q}},t)\,\dot{\mathbf{q}} + \frac{\partial^2 \boldsymbol{\phi}}{\partial t^2}(\mathbf{q},t) = \mathbf{0}. \]

For time-invariant constraints \( \boldsymbol{\phi}(\mathbf{q}) \), the explicit time derivatives vanish, simplifying to

\[ \mathbf{J}_\phi(\mathbf{q})\,\ddot{\mathbf{q}} + \dot{\mathbf{J}}_\phi(\mathbf{q},\dot{\mathbf{q}})\,\dot{\mathbf{q}} = \mathbf{0}. \]

Combining the dynamic and acceleration constraints yields the augmented linear system in unknowns \( (\ddot{\mathbf{q}},\boldsymbol{\lambda}) \):

\[ \begin{bmatrix} \mathbf{M}(\mathbf{q}) & -\mathbf{J}_\phi(\mathbf{q},t)^\top \\ \mathbf{J}_\phi(\mathbf{q},t) & \mathbf{0} \end{bmatrix} \begin{bmatrix} \ddot{\mathbf{q}} \\[0.25em] \boldsymbol{\lambda} \end{bmatrix} = \begin{bmatrix} \mathbf{S}^\top \tau - \mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\,\dot{\mathbf{q}} - \mathbf{g}(\mathbf{q}) \\ - \dot{\mathbf{J}}_\phi(\mathbf{q},\dot{\mathbf{q}},t)\,\dot{\mathbf{q}} - \dfrac{\partial^2 \boldsymbol{\phi}}{\partial t^2}(\mathbf{q},t) \end{bmatrix}. \]

For each configuration \( \mathbf{q} \) and velocity \( \dot{\mathbf{q}} \), solving this linear system yields the constrained accelerations and constraint force multipliers.

As a minimal example, consider a point mass \( m \) moving in the vertical plane with coordinates \( \mathbf{q} = [x\; y]^\top \), constrained to move on a circle of radius \( \ell \) centered at the origin:

\[ \phi(\mathbf{q}) = x^2 + y^2 - \ell^2 = 0. \]

The kinetic and potential energies are

\[ T = \tfrac{1}{2} m (\dot{x}^2 + \dot{y}^2),\quad V = m g y, \]

so the Lagrangian is \( L = T - V \). The constraint Jacobian is

\[ \mathbf{J}_\phi(\mathbf{q}) = \begin{bmatrix} 2x & 2y \end{bmatrix}. \]

The Euler–Lagrange equations with multiplier \( \lambda \) give

\[ \begin{aligned} m \ddot{x} &= 2\lambda x, \\ m \ddot{y} + m g &= 2\lambda y. \end{aligned} \]

The velocity-level constraint is

\[ \mathbf{J}_\phi(\mathbf{q}) \dot{\mathbf{q}} = 2x \dot{x} + 2y \dot{y} = 0, \]

and differentiating once more yields

\[ 2x \ddot{x} + 2y \ddot{y} + 2(\dot{x}^2 + \dot{y}^2) = 0 \quad\Longrightarrow\quad x \ddot{x} + y \ddot{y} + (\dot{x}^2 + \dot{y}^2) = 0. \]

This is precisely the acceleration-level constraint \( \mathbf{J}_\phi \ddot{\mathbf{q}} + \dot{\mathbf{J}}_\phi \dot{\mathbf{q}} = 0 \). Substituting \( \ddot{x},\ddot{y} \) from the dynamic equations into the constraint equation allows solving for \( \lambda \). The resulting \( \lambda \) corresponds to the radial constraint force maintaining constant distance \( \ell \).

This simple example already exhibits the structure that appears in constrained manipulator dynamics: a symmetric inertia matrix, a constraint Jacobian, and a linear algebra problem coupling \( \ddot{\mathbf{q}} \) and \( \boldsymbol{\lambda} \).

For numerical simulation of constrained robot dynamics, a common approach is:

1. Given \( \mathbf{q},\dot{\mathbf{q}} \) at time \( t \), assemble \( \mathbf{M}(\mathbf{q}) \), \( \mathbf{C}(\mathbf{q},\dot{\mathbf{q}}) \), \( \mathbf{g}(\mathbf{q}) \), and \( \mathbf{J}_\phi(\mathbf{q},t) \) (and optionally \( \dot{\mathbf{J}}_\phi \)).
2. Construct the augmented matrix and right-hand side for \( (\ddot{\mathbf{q}},\boldsymbol{\lambda}) \).
3. Solve the linear system, then integrate \( \dot{\mathbf{q}} \) and \( \mathbf{q} \) using a suitable ODE solver.
4. Monitor constraint violations \( \boldsymbol{\phi}(\mathbf{q},t) \) and, if necessary, apply stabilization schemes (discussed in Lesson 4).

The constraint multipliers \( \boldsymbol{\lambda} \) have an important physical interpretation: they are the generalized forces that ensure the generalized coordinates remain on the constraint manifold. In manipulator applications, these are often the joint or contact forces transmitted through closed chains or environmental contacts.

flowchart TD
  S["State (q, qdot, t)"] --> M["Build M(q), C(q,qdot), g(q)"]
  M --> J["Compute J_phi(q,t) and Jdot_phi(q,qdot,t)"]
  J --> AUG["Assemble augmented matrix and rhs"]
  AUG --> SOL["Solve for qddot and lambda"]
  SOL --> INT["Integrate to update (q, qdot)"]
  INT --> CHECK["Check constraint error phi(q,t)"]
  CHECK -->|small| NEXT["Proceed to next step"]
  CHECK -->|large| ADJ["Apply stabilization or smaller time step"]
        

We implement the circle-constrained point mass example using numpy for numerics. Symbolic derivation of \( \mathbf{M} \), \( \mathbf{J}_\phi \), etc., can be automated with sympy or higher-level robotics toolboxes, but here we show the core structure directly.

import numpy as np

def mass_matrix(q, m=1.0):
    # q = [x, y]
    return m * np.eye(2)

def constraint_phi(q, l=1.0):
    x, y = q
    return np.array([x**2 + y**2 - l**2])

def constraint_J(q):
    x, y = q
    return np.array([[2.0 * x, 2.0 * y]])  # shape (1,2)

def constraint_Jdot_qdot(q, qdot):
    x, y = q
    xdot, ydot = qdot
    # J = [2x, 2y]; d/dt J = [2xdot, 2ydot]
    return 2.0 * (xdot**2 + ydot**2)  # scalar Jdot * qdot

def unconstrained_forces(q, qdot, m=1.0, g=9.81):
    # conservative gravity force in generalized coordinates:
    # m yddot = -m g + ...
    # We place gravity as a generalized force on the RHS
    return np.array([0.0, -m * g])

def step_constrained_dynamics(q, qdot, tau=None, m=1.0, l=1.0):
    """
    Compute qddot and lambda for one time instant
    using Lagrange multipliers for a point mass on a circle.
    """
    if tau is None:
        tau = np.zeros(2)

    M = mass_matrix(q, m)
    J = constraint_J(q)
    f = unconstrained_forces(q, qdot, m=m)

    # Build augmented system:
    # [ M  -J^T ] [ qddot ] = [ tau + f ]
    # [ J   0   ] [ lambda]   [  gamma  ]
    # where gamma = -Jdot * qdot for time-invariant constraint
    gamma = -constraint_Jdot_qdot(q, qdot)

    A = np.block([
        [M, -J.T],
        [J, np.zeros((1, 1))]
    ])

    rhs = np.concatenate([tau + f, np.array([gamma])])

    sol = np.linalg.solve(A, rhs)
    qddot = sol[0:2]
    lam = sol[2]

    return qddot, lam

# Example usage:
q = np.array([1.0, 0.0])       # starting at (l, 0)
qdot = np.array([0.0, 1.0])    # tangential velocity
qddot, lam = step_constrained_dynamics(q, qdot)
print("qddot =", qddot)
print("lambda =", lam)
      

In a full simulation, the function step_constrained_dynamics would be called within an ODE right-hand side and integrated using scipy.integrate.solve_ivp or a similar integrator. Robotics-oriented Python libraries (e.g. roboticstoolbox-python) typically provide utilities for computing \( \mathbf{M},\mathbf{C},\mathbf{g},\mathbf{J}_\phi \) for more complex manipulators; the multiplier-based augmentation remains conceptually the same.

In C++, linear algebra for constrained dynamics is commonly handled with Eigen, or with full rigid-body dynamics libraries such as RBDL or Pinocchio which internally use multiplier-like formulations for constraints. Below is a minimal Eigen-based implementation of the same point-mass example:

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

using Eigen::Matrix2d;
using Eigen::Vector2d;
using Eigen::MatrixXd;
using Eigen::VectorXd;

struct CircleSystem {
    double m;
    double l;
    double g;

    CircleSystem(double mass = 1.0, double radius = 1.0, double gravity = 9.81)
        : m(mass), l(radius), g(gravity) {}

    Matrix2d M(const Vector2d& q) const {
        (void)q;
        Matrix2d Mq = Matrix2d::Identity();
        return m * Mq;
    }

    // J_phi(q) = [2x 2y]
    Eigen::RowVector2d J(const Vector2d& q) const {
        return Eigen::RowVector2d(2.0 * q(0), 2.0 * q(1));
    }

    double Jdot_qdot(const Vector2d& q, const Vector2d& qdot) const {
        (void)q;
        double xdot = qdot(0);
        double ydot = qdot(1);
        return 2.0 * (xdot * xdot + ydot * ydot);
    }

    Vector2d unconstrained_forces(const Vector2d& q,
                                  const Vector2d& qdot) const {
        (void)q;
        (void)qdot;
        Vector2d f;
        f << 0.0, -m * g;
        return f;
    }

    // Solve for qddot and lambda
    void step(const Vector2d& q,
              const Vector2d& qdot,
              const Vector2d& tau,
              Vector2d& qddot,
              double& lambda) const
    {
        Matrix2d Mq = M(q);
        Eigen::RowVector2d Jq = J(q);
        double gamma = -Jdot_qdot(q, qdot);
        Vector2d f = unconstrained_forces(q, qdot);

        MatrixXd A(3, 3);
        A.setZero();
        A.block<2,2>(0, 0) = Mq;
        A.block<2,1>(0, 2) = -Jq.transpose();
        A.block<1,2>(2, 0) = Jq;

        VectorXd rhs(3);
        rhs.segment<2>(0) = tau + f;
        rhs(2) = gamma;

        VectorXd sol = A.colPivHouseholderQr().solve(rhs);

        qddot = sol.segment<2>(0);
        lambda = sol(2);
    }
};

int main() {
    CircleSystem sys;
    Vector2d q(1.0, 0.0);
    Vector2d qdot(0.0, 1.0);
    Vector2d tau(0.0, 0.0);
    Vector2d qddot;
    double lambda = 0.0;

    sys.step(q, qdot, tau, qddot, lambda);

    std::cout << "qddot = " << qddot.transpose() << std::endl;
    std::cout << "lambda = " << lambda << std::endl;
    return 0;
}
      

In practical robotics software, the same idea extends to higher-dimensional manipulators by letting M be the manipulator inertia, J the constraint Jacobian (for closed chains or contacts), and solving the augmented system using a robust linear solver.

Modern Java robotics stacks often rely on linear-algebra libraries such as EJML for matrix computations. Below, we illustrate a simple implementation using small fixed-size arrays for the 2D example (sufficient to highlight the multiplier-based formulation). In practice, EJML or similar is recommended for general manipulators.

public class CircleConstrainedDynamics {

    static class State {
        double x, y;
        double xdot, ydot;
    }

    static class Result {
        double xddot, yddot;
        double lambda;
    }

    public static Result step(State s, double m, double g) {
        // M = m * I
        double[][] M = {
            {m, 0.0},
            {0.0, m}
        };

        // J = [2x 2y]
        double[] J = {2.0 * s.x, 2.0 * s.y};

        // Jdot * qdot
        double Jdot_qdot = 2.0 * (s.xdot * s.xdot + s.ydot * s.ydot);
        double gamma = -Jdot_qdot;

        // Unconstrained generalized forces (gravity on y)
        double[] f = {0.0, -m * g};
        double[] tau = {0.0, 0.0};

        // Build augmented 3x3 system A * [xddot, yddot, lambda]^T = rhs
        double[][] A = new double[3][3];

        A[0][0] = M[0][0];
        A[0][1] = M[0][1];
        A[0][2] = -J[0];

        A[1][0] = M[1][0];
        A[1][1] = M[1][1];
        A[1][2] = -J[1];

        A[2][0] = J[0];
        A[2][1] = J[1];
        A[2][2] = 0.0;

        double[] rhs = new double[3];
        rhs[0] = tau[0] + f[0];
        rhs[1] = tau[1] + f[1];
        rhs[2] = gamma;

        // Solve 3x3 linear system via Gaussian elimination (omitted for brevity)
        double[] sol = solve3x3(A, rhs);

        Result r = new Result();
        r.xddot = sol[0];
        r.yddot = sol[1];
        r.lambda = sol[2];
        return r;
    }

    // Simple 3x3 solver (no pivoting)
    static double[] solve3x3(double[][] A, double[] b) {
        double[][] a = new double[3][4];
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) a[i][j] = A[i][j];
            a[i][3] = b[i];
        }
        // Forward elimination
        for (int k = 0; k < 3; ++k) {
            double pivot = a[k][k];
            for (int j = k; j < 4; ++j) a[k][j] /= pivot;
            for (int i = k + 1; i < 3; ++i) {
                double factor = a[i][k];
                for (int j = k; j < 4; ++j) {
                    a[i][j] -= factor * a[k][j];
                }
            }
        }
        // Back substitution
        double[] x = new double[3];
        for (int i = 2; i >= 0; --i) {
            x[i] = a[i][3];
            for (int j = i + 1; j < 3; ++j) {
                x[i] -= a[i][j] * x[j];
            }
        }
        return x;
    }

    public static void main(String[] args) {
        State s = new State();
        s.x = 1.0; s.y = 0.0;
        s.xdot = 0.0; s.ydot = 1.0;

        Result r = step(s, 1.0, 9.81);
        System.out.println("xddot = " + r.xddot + ", yddot = " + r.yddot);
        System.out.println("lambda = " + r.lambda);
    }
}
      

For larger systems, replace the hand-written solver with an EJML-based dense solver, and store \( \mathbf{M} \) and \( \mathbf{J}_\phi \) as EJML matrices.

MATLAB is widely used for robotics dynamics, and the Robotics System Toolbox provides functions for computing \( \mathbf{M},\mathbf{C},\mathbf{g} \) for many manipulator models. Here we show a script-oriented implementation of the circle-constrained example, suitable for embedding in a Simulink block.

function dx = circle_dynamics(t, x, m, l, g)
% x = [q; qdot] = [x; y; xdot; ydot]

q    = x(1:2);
qdot = x(3:4);

M = m * eye(2);
J = [2*q(1), 2*q(2)];

% Jdot * qdot for time-invariant constraint
Jdot_qdot = 2*(qdot(1)^2 + qdot(2)^2);
gamma = -Jdot_qdot;

% Unconstrained generalized forces (gravity)
f = [0; -m*g];
tau = [0; 0];

A = [M, -J'; J, 0];
rhs = [tau + f; gamma];

sol = A \ rhs;
qddot = sol(1:2);  % constrained acceleration

dx = [qdot; qddot];
end

% Example of numerical integration
m = 1.0; l = 1.0; g = 9.81;
x0 = [1.0; 0.0; 0.0; 1.0];  % initial [x; y; xdot; ydot]
[t, x] = ode45(@(t, x) circle_dynamics(t, x, m, l, g), [0 10], x0);
plot(t, x(:,1), t, x(:,2));
legend('x', 'y');
axis equal;
      

In Simulink, the same equations can be implemented using a custom MATLAB Function block, with the augmented matrix solve inside. For general manipulators, one would build \( \mathbf{M},\mathbf{C},\mathbf{g} \) with robotics.RigidBodyTree and add constraint Jacobians for closed chains, then use the multiplier formulation to compute constrained accelerations that feed joint integrator blocks.

Wolfram Mathematica is well suited for symbolic derivation of constrained dynamics. The following notebook-style code derives and simulates the Lagrange-multiplier equations for the circle example:

(* Define symbols *)
Clear["Global`*"];
m = 1; g = 9.81; l = 1;

x[t_]; y[t_];

(* Kinetic and potential energy *)
T = 1/2 m (D[x[t], t]^2 + D[y[t], t]^2);
V = m g y[t];
L = T - V;

(* Holonomic constraint: circle *)
phi = x[t]^2 + y[t]^2 - l^2;
lambda[t_];

(* Euler-Lagrange equations with multiplier *)
eqx = D[D[D[L, D[x[t], t]], t] - D[L, x[t]], t] == 
      lambda[t] * D[phi, x[t]];

eqy = D[D[D[L, D[y[t], t]], t] - D[L, y[t]], t] == 
      lambda[t] * D[phi, y[t]];

(* Constraint equation *)
eqphi = phi == 0;

(* Convert to second-order ODEs *)
vars  = {x[t], y[t], lambda[t]};
dvars = {x''[t], y''[t], lambda'[t]};

(* Solve symbolically for accelerations (could also use DSolve directly) *)
(* For numerical simulation, provide initial conditions and use NDSolve *)
ics = {
  x[0] == l, y[0] == 0,
  x'[0] == 0, y'[0] == 1
};

sol = NDSolve[{eqx, eqy, eqphi, ics}, {x, y, lambda}, {t, 0, 10}];

ParametricPlot[
  Evaluate[{x[t], y[t]} /. sol],
  {t, 0, 10},
  AspectRatio -> 1,
  PlotRange -> All,
  AxesLabel -> {"x", "y"}
]
      

For complex robot mechanisms, symbolic packages can derive \( \mathbf{M},\mathbf{J}_\phi \) and generate optimized code for real-time evaluation, while still relying on the multiplier-based enforcement of constraints.

In this lesson we extended the Lagrange–Euler formulation to handle holonomic constraints via Lagrange multipliers. Starting from d'Alembert’s principle, we showed that ideal constraint forces must lie in the row space of the constraint Jacobian, leading to generalized forces \( \mathbf{J}_\phi^\top\boldsymbol{\lambda} \). For robot manipulators, this yields an augmented system coupling joint accelerations and multipliers, which can be solved at each time step to obtain constrained dynamics. We illustrated the procedure on a simple circle-constrained point mass and implemented the resulting equations in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica. These constructions form the basis for modeling closed-chain mechanisms and contact-constrained robots, which are explored in subsequent lessons.