In robot dynamics, we represent the configuration of an \(n\)-DOF manipulator by a vector of generalized coordinates \( \mathbf{q} = (q_1,\dots,q_n)^\top \). For a serial robot, each \(q_i\) is typically a joint angle (revolute) or displacement (prismatic). The generalized velocities are \( \dot{\mathbf{q}} = (\dot{q}_1,\dots,\dot{q}_n)^\top \).

From Lesson 1, the total kinetic energy of the robot is a quadratic form in the generalized velocities:

\[ T(\mathbf{q},\dot{\mathbf{q}}) = \tfrac{1}{2}\,\dot{\mathbf{q}}^\top \mathbf{M}(\mathbf{q})\,\dot{\mathbf{q}}, \]

where \( \mathbf{M}(\mathbf{q}) \) is the positive-definite joint-space inertia matrix. The total potential energy, due to gravity (and possibly springs), is

\[ V(\mathbf{q}) = V_g(\mathbf{q}) + V_s(\mathbf{q}), \]

where \(V_g\) denotes gravitational potential and \(V_s\) denotes any elastic potential energy (e.g., joint springs, elastic elements).

The Lagrangian is defined as

\[ L(\mathbf{q},\dot{\mathbf{q}},t) = T(\mathbf{q},\dot{\mathbf{q}}) - V(\mathbf{q},t). \]

Note that in most robotic systems, the potential energy may depend on time only through explicitly time-varying elements (e.g., moving base, varying gravity direction), but in this lesson we assume \( V = V(\mathbf{q}) \) for simplicity.

flowchart TD
  A["Robot geometry & inertial data"] --> B["Compute kinetic energy T(q, qd)"]
  A --> C["Compute potential energy V(q)"]
  B --> D["Form L(q, qd) = T - V"]
  C --> D
  D --> E["Apply Euler-Lagrange"]
  E --> F["Equations of motion in q"]
        

The Lagrangian formulation is rooted in a variational principle. Given the Lagrangian \(L(\mathbf{q},\dot{\mathbf{q}},t)\), the action functional over a time interval \([t_0,t_f]\) is

\[ S[\mathbf{q}] = \int_{t_0}^{t_f} L(\mathbf{q}(t),\dot{\mathbf{q}}(t),t)\, \mathrm{d}t. \]

Hamilton's principle states that the actual trajectory \( \mathbf{q}(t) \) followed by the system between fixed endpoints \( \mathbf{q}(t_0) \) and \( \mathbf{q}(t_f) \) makes the action stationary (typically a minimum) with respect to all admissible variations \( \delta\mathbf{q}(t) \) that vanish at the endpoints:

\[ \delta S[\mathbf{q}] = 0 \quad \text{for all variations } \delta\mathbf{q}(t) \text{ with } \delta\mathbf{q}(t_0)=\delta\mathbf{q}(t_f)=\mathbf{0}. \]

This principle generalizes Newton's second law to generalized coordinates and naturally accommodates holonomic constraints that have already been encoded in the choice of \( \mathbf{q} \).

For robotic manipulators, this means:

• We encode rigid-body kinematics into the mapping from joint variables to link positions and orientations.
• From that mapping and inertial parameters, we obtain \(T\) and \(V\).
• We plug \(T\) and \(V\) into the action and impose \( \delta S = 0 \) to derive the joint-space equations of motion.

We now derive the Euler–Lagrange equations for a system with generalized coordinates \( \mathbf{q} = (q_1,\dots,q_n)^\top \). For clarity, we first consider a single coordinate \(q(t)\), then generalize to the vector case.

3.1 Single generalized coordinate

Let \(L(q,\dot{q},t)\) be the Lagrangian of a one-DOF system. Consider a perturbed trajectory \( q_\varepsilon(t) = q(t) + \varepsilon\,\eta(t) \), where \( \eta(t_0) = \eta(t_f) = 0 \) and \( \varepsilon \) is a small scalar. The action for this perturbed trajectory is

\[ S[\varepsilon] = \int_{t_0}^{t_f} L\big(q_\varepsilon(t),\dot{q}_\varepsilon(t),t\big)\,\mathrm{d}t. \]

The stationarity condition is \( \left.\frac{\mathrm{d}S}{\mathrm{d}\varepsilon}\right|\_{\varepsilon=0} = 0 \). Differentiating under the integral sign:

\[ \left.\frac{\mathrm{d}S}{\mathrm{d}\varepsilon}\right|\_{\varepsilon=0} = \int_{t_0}^{t_f} \Bigg( \frac{\partial L}{\partial q}\,\eta(t) + \frac{\partial L}{\partial \dot{q}}\,\dot{\eta}(t) \Bigg)\,\mathrm{d}t. \]

Integrate the second term by parts:

\[ \int_{t_0}^{t_f} \frac{\partial L}{\partial \dot{q}}\,\dot{\eta}(t)\,\mathrm{d}t = \left[ \frac{\partial L}{\partial \dot{q}}\,\eta(t) \right]_{t_0}^{t_f} - \int_{t_0}^{t_f} \frac{\mathrm{d}}{\mathrm{d}t} \left(\frac{\partial L}{\partial \dot{q}}\right)\eta(t)\,\mathrm{d}t. \]

The boundary term vanishes because \( \eta(t_0)=\eta(t_f)=0 \). Thus

\[ \left.\frac{\mathrm{d}S}{\mathrm{d}\varepsilon}\right|\_{\varepsilon=0} = \int_{t_0}^{t_f} \left( \frac{\partial L}{\partial q} - \frac{\mathrm{d}}{\mathrm{d}t} \left(\frac{\partial L}{\partial \dot{q}}\right) \right)\eta(t)\,\mathrm{d}t. \]

Since this must be zero for all admissible \( \eta(t) \), the usual argument from the fundamental lemma of the calculus of variations implies

\[ \frac{\mathrm{d}}{\mathrm{d}t} \left(\frac{\partial L}{\partial \dot{q}}\right) - \frac{\partial L}{\partial q} = 0. \]

3.2 Vector form and generalized forces

For \(n\) generalized coordinates \( \mathbf{q} \), we obtain \(n\) coupled equations:

\[ \frac{\mathrm{d}}{\mathrm{d}t} \left(\frac{\partial L}{\partial \dot{q}_i}\right) - \frac{\partial L}{\partial q_i} = Q_i(\mathbf{q},\dot{\mathbf{q}},t), \quad i=1,\dots,n, \]

where \(Q_i\) are generalized non-conservative forces (from actuators, friction, etc.) introduced via virtual work (as in Lesson 2). In compact vector notation:

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

In future lessons, we will show how these equations can be written in the familiar manipulator form \( \mathbf{M}(\mathbf{q})\ddot{\mathbf{q}} + \mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{q}} + \mathbf{g}(\mathbf{q}) = \boldsymbol{\tau} \), but for now we focus on the Euler–Lagrange structure itself.

4.1 1-DOF mass–spring system

Consider a point mass \(m\) moving along a line with coordinate \(q\) attached to a linear spring of stiffness \(k\) and rest length at \(q=0\). The kinetic and potential energies are

\[ T = \tfrac{1}{2}m\dot{q}^2, \qquad V = \tfrac{1}{2}k q^2. \]

The Lagrangian and its derivatives are

\[ L(q,\dot{q}) = T - V = \tfrac{1}{2}m\dot{q}^2 - \tfrac{1}{2}k q^2, \]

\[ \frac{\partial L}{\partial q} = -k q, \qquad \frac{\partial L}{\partial \dot{q}} = m\dot{q}, \qquad \frac{\mathrm{d}}{\mathrm{d}t} \left(\frac{\partial L}{\partial \dot{q}}\right) = m\ddot{q}. \]

With no non-conservative forces (\(Q=0\)), the Euler–Lagrange equation gives

\[ m\ddot{q} + kq = 0, \]

which is precisely the familiar linear mass–spring equation.

4.2 Planar pendulum (1-DOF) in robotics notation

A classic example closely related to robot joints is the simple planar pendulum with generalized coordinate \(q\) representing the joint angle from the downward vertical. Let the pendulum have length \(l\), mass \(m\) concentrated at the end, and gravity \(g\).

The velocity of the mass is \( v = l\dot{q} \), so the kinetic and potential energies are

\[ T = \tfrac{1}{2} m (l\dot{q})^2 = \tfrac{1}{2} m l^2 \dot{q}^2, \qquad V = m g l (1 - \cos q). \]

(We choose \(V=0\) at the bottom position.) Then

\[ L(q,\dot{q}) = \tfrac{1}{2} m l^2 \dot{q}^2 - m g l (1 - \cos q). \]

Derivatives:

\[ \frac{\partial L}{\partial q} = -m g l \sin q, \qquad \frac{\partial L}{\partial \dot{q}} = m l^2 \dot{q}, \qquad \frac{\mathrm{d}}{\mathrm{d}t} \left(\frac{\partial L}{\partial \dot{q}}\right) = m l^2 \ddot{q}. \]

With generalized torque \(Q = \tau\) applied at the joint, the Euler–Lagrange equation yields

\[ m l^2 \ddot{q} + m g l \sin q = \tau. \]

This is the standard equation of motion of a 1-DOF robot joint under gravity (without friction). In further lessons we will generalize this to multiple joints and identify the inertia, Coriolis/centrifugal, and gravity terms systematically.

For a serial robot with known link parameters, the Lagrangian-based modeling procedure can be summarized as:

1. Choose generalized coordinates \( \mathbf{q} \) (typically joint variables).
2. Use forward kinematics (from previous chapters) to compute each link center-of-mass position and orientation as functions of \( \mathbf{q} \).
3. Compute linear and angular velocities of each link and assemble the total kinetic energy \(T\).
4. Compute gravitational and elastic potential energies to obtain \(V\).
5. Form the Lagrangian \(L = T - V\).
6. Compute \( \frac{\partial L}{\partial \mathbf{q}} \) and \( \frac{\partial L}{\partial \dot{\mathbf{q}}} \), then \( \frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\partial L}{\partial \dot{\mathbf{q}}}\right) \).
7. Include generalized input torques/forces and non-conservative forces via \( \mathbf{Q} \), and write the Euler–Lagrange equations.

flowchart TD
  S["Start: choose joint coordinates q"] --> FK["Kinematics: positions & velocities of links"]
  FK --> TE["Compute kinetic energy T(q, qd)"]
  FK --> PE["Compute potential energy V(q)"]
  TE --> L["Form L(q, qd) = T - V"]
  PE --> L
  L --> EL["Compute d/dt(dL/dqd) - dL/dq"]
  EL --> QF["Add generalized forces Q(q, qd, t)"]
  QF --> EOM["Final equations of motion"]
        

In Python, symbolic derivation of Lagrange equations is conveniently handled with sympy (and optionally sympy.physics.mechanics). For robotics, high-level libraries such as roboticstoolbox-python build on these ideas but often start from pre-derived dynamic models.

Below we derive the planar pendulum equation symbolically and then build a numeric right-hand side suitable for simulation.

import sympy as sp

# Symbols
t = sp.symbols('t')
q = sp.Function('q')(t)          # generalized coordinate q(t)
m, l, g, tau = sp.symbols('m l g tau')

# Kinetic and potential energies
qd = sp.diff(q, t)
T = sp.Rational(1, 2) * m * l**2 * qd**2
V = m * g * l * (1 - sp.cos(q))

L = T - V

# Euler-Lagrange equation: d/dt(dL/dqd) - dL/dq = tau
dL_dq  = sp.diff(L, q)
dL_dqd = sp.diff(L, qd)
d_dt_dL_dqd = sp.diff(dL_dqd, t)

EL_eq = sp.simplify(d_dt_dL_dqd - dL_dq - tau)

print("Euler-Lagrange equation (symbolic):")
print(EL_eq)  # should be m*l**2*qdd + m*g*l*sp.sin(q) - tau

# Replace q''(t) symbolically for clarity
qdd = sp.Function('qdd')(t)
EL_eq_sub = EL_eq.subs(sp.diff(q, (t, 2)), qdd)
print("With qdd substitution:")
print(EL_eq_sub)

# Build a numeric function for simulation using lambdify
q_sym, qd_sym, qdd_sym = sp.symbols('q_sym qd_sym qdd_sym')
EL_numeric = m*l**2 * qdd_sym + m*g*l*sp.sin(q_sym) - tau

rhs_qdd = sp.solve(sp.Eq(EL_numeric, 0), qdd_sym)[0]

# Create a Python function for qdd(q, qd, tau)
qdd_fun = sp.lambdify((q_sym, qd_sym, tau, m, l, g), rhs_qdd, 'numpy')

# Example usage inside an ODE integrator (e.g., scipy.integrate.solve_ivp):
def pendulum_rhs(t, x, m_val, l_val, g_val, tau_fun):
    q_val, qd_val = x
    tau_val = tau_fun(t, q_val, qd_val)  # user-defined torque policy
    qdd_val = qdd_fun(q_val, qd_val, tau_val, m_val, l_val, g_val)
    return [qd_val, qdd_val]
      

This approach generalizes to multi-DOF robots by using vector coordinates \( \mathbf{q} \) and constructing the total kinetic and potential energies from link kinematics and inertias.

In C++, Lagrangian dynamics for robots are commonly handled via libraries such as RBDL, Pinocchio, or Orocos KDL, often combined with Eigen for linear algebra. Here we implement the 1-DOF pendulum equation \( m l^2 \ddot{q} + m g l \sin q = \tau \) directly.

#include <cmath>
#include <array>

// Simple pendulum dynamics: x = [q, qd]
struct PendulumParams {
  double m;
  double l;
  double g;
};

// User-specified torque profile
double torque(double t, double q, double qd) {
  // Example: simple PD control about q = 0
  double Kp = 5.0;
  double Kd = 1.0;
  return -Kp * q - Kd * qd;
}

// Right-hand side of the ODE: xdot = f(t, x)
std::array<double, 2> pendulum_rhs(double t,
                                    const std::array<double, 2>& x,
                                    const PendulumParams& p) {
  double q  = x[0];
  double qd = x[1];

  double tau = torque(t, q, qd);

  double denom = p.m * p.l * p.l;
  double qdd   = (tau - p.m * p.g * p.l * std::sin(q)) / denom;

  return {qd, qdd};
}
      

For a full robot model, libraries like RBDL or Pinocchio internally construct the Lagrangian-based dynamics from link descriptions (URDF) and provide functions that directly compute \( \mathbf{M}(\mathbf{q}) \), \( \mathbf{C}(\mathbf{q},\dot{\mathbf{q}}) \), and \( \mathbf{g}(\mathbf{q}) \) rather than implementing the Euler–Lagrange derivation manually.

Java does not have a dominant robotics dynamics library, but we can implement Lagrangian equations using general-purpose matrix libraries such as EJML or Apache Commons Math. The following example shows the 1-DOF pendulum dynamics using plain Java.

public class Pendulum {

    public static class Params {
        public double m;
        public double l;
        public double g;

        public Params(double m, double l, double g) {
            this.m = m;
            this.l = l;
            this.g = g;
        }
    }

    // Example torque profile (PD control)
    public static double torque(double t, double q, double qd) {
        double Kp = 5.0;
        double Kd = 1.0;
        return -Kp * q - Kd * qd;
    }

    // x = [q, qd], returns [qd, qdd]
    public static double[] rhs(double t, double[] x, Params p) {
        double q  = x[0];
        double qd = x[1];

        double tau = torque(t, q, qd);

        double denom = p.m * p.l * p.l;
        double qdd   = (tau - p.m * p.g * p.l * Math.sin(q)) / denom;

        return new double[]{qd, qdd};
    }
}
      

For multi-DOF robots, you would represent vectors \( \mathbf{q} \), \( \dot{\mathbf{q}} \), and matrices \( \mathbf{M}(\mathbf{q}) \), etc., using a matrix library, and implement the Lagrangian-based equations in joint space exactly as in the scalar example but in vector form.

MATLAB, together with the Robotics System Toolbox and Simscape Multibody, is widely used for Lagrangian robot dynamics. Here we implement the pendulum equation in MATLAB code and outline a basic Simulink model.

function dx = pendulum_rhs(t, x, params)
% x = [q; qd]
q  = x(1);
qd = x(2);

m = params.m;
l = params.l;
g = params.g;

% Example torque: PD about q = 0
Kp = 5.0;
Kd = 1.0;
tau = -Kp * q - Kd * qd;

qdd = (tau - m * g * l * sin(q)) / (m * l^2);

dx = [qd; qdd];
end

% Example of numerical integration:
params.m = 1.0;
params.l = 1.0;
params.g = 9.81;

x0 = [0.2; 0.0];  % initial angle and velocity
tspan = [0 10];

[t_sol, x_sol] = ode45(@(t, x) pendulum_rhs(t, x, params), tspan, x0);
      

In Simulink, you can:

• Create an integrator chain for \(q\) and \(\dot{q}\).
• Compute \( \tau \) from \(q\) and \(\dot{q}\) (e.g., via a Subsystem implementing a PD control law).
• Implement \( \ddot{q} = ( \tau - m g l \sin(q) ) / (m l^2) \) in a block diagram.
• Use scopes to visualize \(q(t)\) and \(\dot{q}(t)\).

Wolfram Mathematica is well suited to symbolic derivation of Lagrangian dynamics, with built-in differentiation and equation-solving tools.

(* Define symbols and generalized coordinate *)
ClearAll["Global`*"];
m  = Symbol["m"];
l  = Symbol["l"];
g  = Symbol["g"];
tau = Symbol["tau"];

q  = q[t];
qd = D[q, t];

(* Kinetic and potential energies *)
T = 1/2 m l^2 qd^2;
V = m g l (1 - Cos[q]);

L = T - V;

(* Euler-Lagrange equation *)
eq = D[D[L, qd], t] - D[L, q] == tau;

(* Simplify and solve for q''(t) *)
eqSimplified = Simplify[eq];
qdd = D[q, {t, 2}];

sol = Solve[eqSimplified, qdd][[1, 1]];

Print["Equation of motion:"];
Print[eqSimplified];
Print["qdd(t) = ", qdd /. sol];
      

For multi-DOF manipulators, you would introduce a list of generalized coordinates \(\{q_1(t),\dots,q_n(t)\}\), compute total \(T\) and \(V\), and then generate \(n\) Euler–Lagrange equations systematically with a loop over coordinates.

In this lesson we introduced the Lagrangian formulation of robot dynamics. Starting from generalized coordinates and the kinetic and potential energy expressions of multi-body systems, we constructed the Lagrangian \(L = T - V\) and applied Hamilton's principle to derive the Euler–Lagrange equations \( \frac{\mathrm{d}}{\mathrm{d}t} \left(\frac{\partial L}{\partial \dot{\mathbf{q}}}\right) - \frac{\partial L}{\partial \mathbf{q}} = \mathbf{Q} \). We illustrated these equations on standard robotic examples (mass–spring, pendulum, coupled masses), emphasized the coordinate-invariant nature of the formulation, and implemented basic Lagrangian dynamics in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica. In the next lesson, we will structure these equations into the standard manipulator form and analyze the properties of the inertia, Coriolis, and gravity terms.