Recall from the previous lessons that for an \(n\)-DOF serial manipulator, the Lagrange–Euler equations can be written in the standard form

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

where \( \mathbf{q} \in \mathbb{R}^n \) are generalized joint coordinates, \( \mathbf{M}(\mathbf{q}) \) is the inertia matrix, \( \mathbf{C}(\mathbf{q},\dot{\mathbf{q}}) \) groups Coriolis and centrifugal effects, \( \mathbf{g}(\mathbf{q}) \) is the gravity vector, and \( \boldsymbol{\tau} \) (written in MathJax as \(\tau\)) is the vector of joint torques/forces.

In this lesson we will:

• Define kinetic, potential, and total mechanical energy of a manipulator.
• Show that the Lagrange–Euler equations imply an exact power balance: \( \dot{H} = \dot{\mathbf{q}}^\top \boldsymbol{\tau} \) for conservative systems.
• Derive the skew-symmetry property of \( \dot{\mathbf{M}}(\mathbf{q}) - 2\mathbf{C}(\mathbf{q},\dot{\mathbf{q}}) \) and relate it to energy consistency.
• Introduce passivity and show that robot manipulators are passive from input \( \boldsymbol{\tau} \) to output \( \dot{\mathbf{q}} \).
• Develop numerical procedures and code to test these properties for a concrete model.

flowchart TD
  TQ["Joint torques tau(t)"] --> DYN["Dynamics: M(q) qddot + C(q,qd) qd + g(q) = tau"]
  DYN --> VEL["Joint velocities qd(t)"]
  DYN --> ACC["Joint accelerations qddot(t)"]
  VEL --> POS["Joint positions q(t)"]
  VEL --> KE["Kinetic energy T(q,qd)"]
  POS --> PE["Potential energy U(q)"]
  KE --> H["Total energy H = T + U"]
  PE --> H
  H --> PWR["Check dH/dt vs power tau^T qd"]
        

The kinetic energy of an \(n\)-DOF manipulator is quadratic in joint 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}) \) arises from summing link translational and rotational energies. Fundamental properties (without proof here, but following from rigid-body mechanics) are:

\[ \mathbf{M}(\mathbf{q}) = \mathbf{M}(\mathbf{q})^\top,\qquad \mathbf{v}^\top \mathbf{M}(\mathbf{q}) \mathbf{v} > 0,\ \forall \mathbf{v} \neq \mathbf{0}. \]

Thus \( \mathbf{M}(\mathbf{q}) \) is symmetric positive definite for all physically valid \( \mathbf{q} \).

The potential energy typically consists of gravitational and possibly elastic contributions:

\[ U(\mathbf{q}) = U_g(\mathbf{q}) + U_{\text{elastic}}(\mathbf{q}), \]

and the total mechanical energy is

\[ H(\mathbf{q},\dot{\mathbf{q}}) = T(\mathbf{q},\dot{\mathbf{q}}) + U(\mathbf{q}). \]

In Lagrangian mechanics we define the Lagrangian \( L(\mathbf{q},\dot{\mathbf{q}}) = T(\mathbf{q},\dot{\mathbf{q}}) - U(\mathbf{q}) \). Gravity enters the equations via \( \mathbf{g}(\mathbf{q}) = \nabla_{\mathbf{q}} U(\mathbf{q}) \), consistent with the generalized forces generated by potential fields.

Starting from the standard dynamics

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

multiply on the left by \( \dot{\mathbf{q}}^\top \):

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

On the other hand, differentiate the kinetic energy:

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

Subtracting these two expressions eliminates the term \( \dot{\mathbf{q}}^\top \mathbf{M}(\mathbf{q}) \ddot{\mathbf{q}} \) and yields

\[ \dot{T} - \dot{\mathbf{q}}^\top \boldsymbol{\tau} + \dot{\mathbf{q}}^\top \mathbf{g}(\mathbf{q}) = \tfrac{1}{2}\dot{\mathbf{q}}^\top \dot{\mathbf{M}}(\mathbf{q}) \dot{\mathbf{q}} - \dot{\mathbf{q}}^\top \mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{q}}. \]

This shows that power balance will hold if the right-hand side vanishes; i.e., if

\[ \dot{\mathbf{q}}^\top \left( \tfrac{1}{2}\dot{\mathbf{M}}(\mathbf{q}) - \mathbf{C}(\mathbf{q},\dot{\mathbf{q}}) \right) \dot{\mathbf{q}} = 0 \quad \forall\ \mathbf{q},\dot{\mathbf{q}}. \]

Define the matrix

\[ \mathbf{S}(\mathbf{q},\dot{\mathbf{q}}) := \dot{\mathbf{M}}(\mathbf{q}) - 2\mathbf{C}(\mathbf{q},\dot{\mathbf{q}}). \]

A sufficient and necessary condition for \( \dot{\mathbf{q}}^\top \left(\tfrac{1}{2}\dot{\mathbf{M}} - \mathbf{C}\right)\dot{\mathbf{q}} = 0 \) for all \( \dot{\mathbf{q}} \) is that \( \mathbf{S} \) be skew-symmetric:

\[ \mathbf{S}(\mathbf{q},\dot{\mathbf{q}}) = -\mathbf{S}(\mathbf{q},\dot{\mathbf{q}})^\top \quad\Longleftrightarrow\quad \mathbf{x}^\top \mathbf{S}(\mathbf{q},\dot{\mathbf{q}})\mathbf{x} = 0 \ \forall \mathbf{x}. \]

Under this condition,

\[ \dot{T} = \dot{\mathbf{q}}^\top \boldsymbol{\tau} - \dot{\mathbf{q}}^\top \mathbf{g}(\mathbf{q}). \]

Since \( \mathbf{g}(\mathbf{q}) = \nabla_{\mathbf{q}} U(\mathbf{q}) \), the time derivative of the potential energy is

\[ \dot{U} = \nabla_{\mathbf{q}} U(\mathbf{q})^\top \dot{\mathbf{q}} = \mathbf{g}(\mathbf{q})^\top \dot{\mathbf{q}}. \]

Therefore, the total energy derivative is

\[ \dot{H} = \dot{T} + \dot{U} = \dot{\mathbf{q}}^\top \boldsymbol{\tau}, \]

which states that the rate of change of mechanical energy equals the instantaneous mechanical power injected by the joint torques. This identity is the core energy consistency check for any candidate \( \mathbf{M},\mathbf{C},\mathbf{g} \).

Consider a system with state \( \mathbf{x} \), input \( \mathbf{u} \), and output \( \mathbf{y} \). A nonnegative function \( S(\mathbf{x}) \) is a storage function if it upper-bounds the internal energy. The system is passive from \( \mathbf{u} \) to \( \mathbf{y} \) if

\[ S(\mathbf{x}(t)) - S(\mathbf{x}(0)) \le \int_0^t \mathbf{y}(\sigma)^\top \mathbf{u}(\sigma)\,\mathrm{d}\sigma \quad \forall t \ge 0, \]

for all trajectories and all inputs. If equality holds for all trajectories, the system is lossless passive. If

\[ S(\mathbf{x}(t)) - S(\mathbf{x}(0)) \le \int_0^t \mathbf{y}(\sigma)^\top \mathbf{u}(\sigma)\,\mathrm{d}\sigma - \alpha \int_0^t \|\mathbf{y}(\sigma)\|^2\,\mathrm{d}\sigma \quad (\alpha > 0), \]

the system is output strictly passive.

For a robot manipulator, a natural choice of port variables is \( \mathbf{u} = \boldsymbol{\tau} \) and \( \mathbf{y} = \dot{\mathbf{q}} \), with storage \( S = H \) equal to the total mechanical energy. If the energy balance \( \dot{H} = \dot{\mathbf{q}}^\top \boldsymbol{\tau} \) holds, then integrating in time shows that the manipulator is lossless passive.

Including viscous joint friction with positive-semidefinite damping matrix \( \mathbf{D} \), the dynamics become

\[ \mathbf{M}(\mathbf{q}) \ddot{\mathbf{q}} + \mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{q}} + \mathbf{g}(\mathbf{q}) + \mathbf{D}\dot{\mathbf{q}} = \boldsymbol{\tau}. \]

Repeating the previous derivation yields

\[ \dot{H} = \dot{\mathbf{q}}^\top \boldsymbol{\tau} - \dot{\mathbf{q}}^\top \mathbf{D} \dot{\mathbf{q}}. \]

Since \( \mathbf{D} \) is positive semidefinite, \( \dot{\mathbf{q}}^\top \mathbf{D}\dot{\mathbf{q}} \ge 0 \), so

\[ \dot{H} \le \dot{\mathbf{q}}^\top \boldsymbol{\tau}. \]

Hence the manipulator with viscous damping is output strictly passive from \( \boldsymbol{\tau} \) to \( \dot{\mathbf{q}} \). This crucial property is often exploited in energy-based analyses and, later, in control design (although we will not introduce control laws in this chapter).

flowchart TD
  U["Input: joint torques tau"] --> SYS["Mechanical dynamics"]
  SYS --> Y["Output: joint velocities qd"]
  SYS --> HNODE["Energy H = T + U"]
  HNODE --> DISS["Damping: qd^T D qd >= 0"]
  U --> PINT["Power tau^T qd"]
  PINT --> HNODE
  DISS --> HDEC["Energy non-increasing for tau = 0"]
        

Suppose you have derived symbolic or numeric expressions for \( \mathbf{M}(\mathbf{q}) \), \( \mathbf{C}(\mathbf{q},\dot{\mathbf{q}}) \), and \( \mathbf{g}(\mathbf{q}) \). A practical procedure to test energy consistency is:

1. Check symmetry and positive definiteness of \( \mathbf{M}(\mathbf{q}) \).
   For random joint configurations \( \mathbf{q} \), numerically compute eigenvalues of \( \mathbf{M}(\mathbf{q}) \). All eigenvalues must be strictly positive and \( \|\mathbf{M} - \mathbf{M}^\top\| \) must be near numerical zero.
2. Check skew-symmetry of \( \dot{\mathbf{M}} - 2\mathbf{C} \).
   For random \( \mathbf{q},\dot{\mathbf{q}} \), approximate the directional derivative of \( \mathbf{M} \) along \( \dot{\mathbf{q}} \):

   \[ \dot{\mathbf{M}}(\mathbf{q})\dot{\mathbf{q}} \approx \frac{\mathbf{M}(\mathbf{q} + \varepsilon \dot{\mathbf{q}}) - \mathbf{M}(\mathbf{q} - \varepsilon \dot{\mathbf{q}})} {2\varepsilon} \dot{\mathbf{q}}, \]

   then evaluate the scalar

   \[ \sigma = \dot{\mathbf{q}}^\top \left( \tfrac{1}{2}\dot{\mathbf{M}}(\mathbf{q}) - \mathbf{C}(\mathbf{q},\dot{\mathbf{q}}) \right) \dot{\mathbf{q}}. \]

   For a consistent model, \( \sigma \) should be numerically close to zero.
3. Check power balance along simulated trajectories.
   Choose an input \( \boldsymbol{\tau}(t) \), simulate the dynamics, and numerically compute:

   \[ H(t_k) = \tfrac{1}{2}\dot{\mathbf{q}}(t_k)^\top \mathbf{M}(\mathbf{q}(t_k)) \dot{\mathbf{q}}(t_k) + U(\mathbf{q}(t_k)), \]

   together with an approximate derivative \( \dot{H}(t_k) \) via finite differences, and compare to the power

   \[ P(t_k) = \dot{\mathbf{q}}(t_k)^\top \boldsymbol{\tau}(t_k). \]

   For a frictionless model, \( \dot{H}(t_k) \approx P(t_k) \); with damping, \( \dot{H}(t_k) \le P(t_k) \).
4. Zero-input consistency.
   With \( \boldsymbol{\tau}(t) \equiv \mathbf{0} \) and no damping, total energy should be numerically conserved; with damping, it should be non-increasing.

flowchart TD
  A["Derived M(q), C(q,qd), g(q)"] --> B["Check M symmetric / SPD"]
  B --> C["Test v^T (dM_dt - 2 C) v approx 0"]
  C --> D["Simulate trajectories with tau(t)"]
  D --> E["Compute H(t) and dH/dt"]
  E --> F["Compare dH/dt with tau^T qd"]
  F --> G["Model passes energy and passivity checks?"]
        

In Python, we can use numpy for linear algebra and either implement \( \mathbf{M}, \mathbf{C}, \mathbf{g}, U \) by hand or call robotics libraries such as roboticstoolbox-python or pinocchio (not shown here) to obtain them. Below is a template that assumes you already have functions M(q), C(q, dq), g(q), and U(q) for your robot.

import numpy as np

# User-supplied robot model -----------------------------------------

def M(q):
    """
    Inertia matrix M(q) for your manipulator.
    q: ndarray shape (n,)
    Returns ndarray shape (n, n)
    """
    raise NotImplementedError("Implement M(q) for your robot")

def C(q, dq):
    """
    Coriolis/centrifugal matrix C(q, dq).
    """
    raise NotImplementedError("Implement C(q, dq) for your robot")

def g_vec(q):
    """
    Gravity vector g(q) = grad_q U(q).
    """
    raise NotImplementedError("Implement g(q) for your robot")

def potential_energy(q):
    """
    Potential energy U(q).
    """
    raise NotImplementedError("Implement U(q) for your robot")

# Energy-related functions ------------------------------------------

def kinetic_energy(q, dq):
    Mq = M(q)
    return 0.5 * dq.T @ Mq @ dq

def total_energy(q, dq):
    return kinetic_energy(q, dq) + potential_energy(q)

def directional_M_dot(q, dq, eps=1e-6):
    """
    Approximate dot(M) along direction dq using a central difference.
    Returns a matrix of the same shape as M(q).
    """
    q_plus = q + eps * dq
    q_minus = q - eps * dq
    M_plus = M(q_plus)
    M_minus = M(q_minus)
    # Directional derivative: d/dt M(q(t)) with qdot = dq
    return (M_plus - M_minus) / (2.0 * eps)

def scalar_sigma(q, dq):
    """
    Compute sigma = dq^T (0.5 * dot(M) - C) dq, which should be ~0.
    """
    Mdot = directional_M_dot(q, dq)
    Cmat = C(q, dq)
    middle = 0.5 * Mdot - Cmat
    return dq.T @ (middle @ dq)

# Simple time stepping to test energy balance -----------------------

def dynamics_rhs(q, dq, tau):
    """
    Compute qdd given q, dq, tau using M(q) qdd + C(q,dq) dq + g(q) = tau.
    """
    Mq = M(q)
    Cq = C(q, dq)
    gq = g_vec(q)
    rhs = tau - Cq @ dq - gq
    # Solve M qdd = rhs
    qdd = np.linalg.solve(Mq, rhs)
    return qdd

def simulate_step(q, dq, tau, dt):
    """
    One step of explicit Euler integration for testing only.
    """
    qdd = dynamics_rhs(q, dq, tau)
    q_next = q + dt * dq
    dq_next = dq + dt * qdd
    return q_next, dq_next

def energy_power_check(q0, dq0, tau_func, dt=1e-3, steps=1000):
    """
    Simulate and record H(t) and power dq^T tau(t).
    tau_func(t, q, dq) should return ndarray of shape (n,).
    """
    q = q0.copy()
    dq = dq0.copy()
    H_vals = []
    P_vals = []
    t_vals = []

    for k in range(steps):
        t = k * dt
        tau = tau_func(t, q, dq)
        H_vals.append(total_energy(q, dq))
        P_vals.append(dq.T @ tau)
        t_vals.append(t)
        q, dq = simulate_step(q, dq, tau, dt)

    return np.array(t_vals), np.array(H_vals), np.array(P_vals)

# Example usage (after implementing M, C, g, U):
# q0 = np.zeros(n)
# dq0 = np.zeros(n)
# def zero_tau(t, q, dq):
#     return np.zeros_like(q)
# t, H, P = energy_power_check(q0, dq0, zero_tau)
# print("Energy variation:", H.max() - H.min())
# print("Mean power:", P.mean())
      

This code allows you to numerically verify:

• Near-zero values of scalar_sigma(q, dq) for random states.
• Approximate conservation of total_energy for zero input and no damping.
• Consistency between the time derivative of total_energy and the power dq.T @ tau for nonzero inputs.

In C++, a common approach is to represent vectors and matrices using the Eigen library, and to use dynamics libraries (e.g. RBDL, Pinocchio C++ bindings) or your own implementation for \( \mathbf{M}, \mathbf{C}, \mathbf{g} \). Below is a minimal skeleton focusing on energy checks:

#include <Eigen/Dense>
#include <functional>

// Type aliases
using Vector = Eigen::VectorXd;
using Matrix = Eigen::MatrixXd;

// User-supplied dynamics (to be implemented or connected to a robotics library)
Matrix M(const Vector& q);
Matrix C(const Vector& q, const Vector& dq);
Vector g_vec(const Vector& q);
double potentialEnergy(const Vector& q);

double kineticEnergy(const Vector& q, const Vector& dq) {
    Matrix Mq = M(q);
    return 0.5 * dq.transpose() * Mq * dq;
}

double totalEnergy(const Vector& q, const Vector& dq) {
    return kineticEnergy(q, dq) + potentialEnergy(q);
}

Matrix directionalMDot(const Vector& q, const Vector& dq, double eps = 1e-6) {
    Vector q_plus = q + eps * dq;
    Vector q_minus = q - eps * dq;
    Matrix M_plus = M(q_plus);
    Matrix M_minus = M(q_minus);
    return (M_plus - M_minus) / (2.0 * eps);
}

double scalarSigma(const Vector& q, const Vector& dq) {
    Matrix Mdot = directionalMDot(q, dq);
    Matrix Cmat = C(q, dq);
    Matrix middle = 0.5 * Mdot - Cmat;
    return dq.transpose() * middle * dq;
}

Vector dynamicsRhs(const Vector& q, const Vector& dq, const Vector& tau) {
    Matrix Mq = M(q);
    Matrix Cq = C(q, dq);
    Vector gq = g_vec(q);
    Vector rhs = tau - Cq * dq - gq;
    // Solve M qdd = rhs
    return Mq.ldlt().solve(rhs);
}

void stepEuler(Vector& q, Vector& dq,
               const std::function<Vector(double, const Vector&, const Vector&)>& tauFunc,
               double t, double dt) {
    Vector tau = tauFunc(t, q, dq);
    Vector qdd = dynamicsRhs(q, dq, tau);
    q += dt * dq;
    dq += dt * qdd;
}
      

Using this skeleton, you can run Monte Carlo tests of scalarSigma(q, dq), energy conservation for tau = 0, or power balance for arbitrary torque trajectories.

Java does not have a de facto standard robotics dynamics library, but linear algebra libraries like EJML can be used. For clarity, we show a plain double-array implementation of the key checks, assuming the user provides M(q), C(q,dq), g(q), and U(q) in compatible formats.

public class EnergyChecks {

    // User-supplied functions (placeholders)
    // M(q): returns n x n matrix
    public static double[][] M(double[] q) {
        throw new UnsupportedOperationException("Implement M(q)");
    }

    // C(q,dq): returns n x n matrix
    public static double[][] C(double[] q, double[] dq) {
        throw new UnsupportedOperationException("Implement C(q,dq)");
    }

    // g(q): returns n-vector
    public static double[] g(double[] q) {
        throw new UnsupportedOperationException("Implement g(q)");
    }

    public static double potentialEnergy(double[] q) {
        throw new UnsupportedOperationException("Implement U(q)");
    }

    // Basic linear algebra helpers --------------------------------------

    public static double[] matVec(double[][] A, double[] x) {
        int n = x.length;
        double[] y = new double[n];
        for (int i = 0; i < n; ++i) {
            double sum = 0.0;
            for (int j = 0; j < n; ++j) {
                sum += A[i][j] * x[j];
            }
            y[i] = sum;
        }
        return y;
    }

    public static double dot(double[] a, double[] b) {
        double s = 0.0;
        for (int i = 0; i < a.length; ++i) {
            s += a[i] * b[i];
        }
        return s;
    }

    public static double kineticEnergy(double[] q, double[] dq) {
        double[][] Mq = M(q);
        double[] tmp = matVec(Mq, dq);
        return 0.5 * dot(dq, tmp);
    }

    public static double totalEnergy(double[] q, double[] dq) {
        return kineticEnergy(q, dq) + potentialEnergy(q);
    }

    public static double[][] directionalMDot(double[] q, double[] dq, double eps) {
        int n = q.length;
        double[] qPlus = new double[n];
        double[] qMinus = new double[n];
        for (int i = 0; i < n; ++i) {
            qPlus[i] = q[i] + eps * dq[i];
            qMinus[i] = q[i] - eps * dq[i];
        }
        double[][] Mplus = M(qPlus);
        double[][] Mminus = M(qMinus);
        double[][] Mdot = new double[n][n];
        double inv = 1.0 / (2.0 * eps);
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) {
                Mdot[i][j] = (Mplus[i][j] - Mminus[i][j]) * inv;
            }
        }
        return Mdot;
    }

    public static double scalarSigma(double[] q, double[] dq) {
        int n = q.length;
        double[][] Mdot = directionalMDot(q, dq, 1e-6);
        double[][] Cmat = C(q, dq);
        double[][] middle = new double[n][n];
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) {
                middle[i][j] = 0.5 * Mdot[i][j] - Cmat[i][j];
            }
        }
        double[] tmp = matVec(middle, dq);
        return dot(dq, tmp);
    }
}
      

This structure can be embedded into a larger simulation loop to compare energy variation and power \( \dot{\mathbf{q}}^\top \boldsymbol{\tau} \) over time.

In MATLAB, robot dynamics are often handled by toolboxes such as Robotics System Toolbox or Peter Corke's Robotics Toolbox. Here we show generic functions for energy and skew-symmetry checks. They can be called from a script or from a Simulink block via MATLAB Function blocks.

function H = total_energy(q, dq)
    T = kinetic_energy(q, dq);
    U = potential_energy(q);
    H = T + U;
end

function T = kinetic_energy(q, dq)
    Mq = M_matrix(q);
    T = 0.5 * dq.' * Mq * dq;
end

function sigma = scalar_sigma(q, dq)
    Mdot = directional_M_dot(q, dq);
    Cq = C_matrix(q, dq);
    middle = 0.5 * Mdot - Cq;
    sigma = dq.' * middle * dq;
end

function Mdot = directional_M_dot(q, dq)
    eps = 1e-6;
    q_plus = q + eps * dq;
    q_minus = q - eps * dq;
    M_plus = M_matrix(q_plus);
    M_minus = M_matrix(q_minus);
    Mdot = (M_plus - M_minus) / (2 * eps);
end

% User must implement:
%   M_matrix(q)  - inertia matrix
%   C_matrix(q,dq) - Coriolis matrix
%   g_vector(q) - gravity vector
%   potential_energy(q) - potential
      

In Simulink, one can build a subsystem that:

• Receives q, dq, and tau as signals.
• Computes H(t) using total_energy(q,dq).
• Computes instantaneous power dq.'*tau.
• Scopes both signals to visually verify \( \dot{H}(t) \approx \dot{\mathbf{q}}^\top \boldsymbol{\tau} \).

Mathematica is well-suited for symbolic derivations and exact verification of skew-symmetry and passivity properties for low-DOF manipulators.

(* q is a vector of generalized coordinates, dq a vector of velocities *)
(* Define symbolic M(q), C(q,dq), g(q) and potential U(q) *)

(* Example: treat M as a matrix-valued function of q *)
M[q_] := { {m11[q], m12[q]}, {m12[q], m22[q]} }; (* 2-DOF example structure *)
C[q_, dq_] := { {c11[q, dq], c12[q, dq]}, {c21[q, dq], c22[q, dq]} };
gVec[q_] := {g1[q], g2[q]};
U[q_] := ug[q];  (* scalar potential *)

kineticEnergy[q_, dq_] := 1/2 dq . (M[q].dq);
totalEnergy[q_, dq_] := kineticEnergy[q, dq] + U[q];

(* Time derivatives via chain rule: q[t], dq[t] = q'[t] *)
Clear[q1, q2];
q1[t_]; q2[t_];
qSym[t_] := {q1[t], q2[t]};
dqSym[t_] := D[qSym[t], t];

(* Compute d/dt of kinetic energy along trajectories *)
Tsym[t_] := kineticEnergy[qSym[t], dqSym[t]];
dTdt = D[Tsym[t], t] // Simplify;

(* Define S(q, dq) = dot(M) - 2 C and check skew-symmetry condition x^T S x == 0 *)
(* For exact check, expand x^T S x symbolically and simplify *)

x = {x1, x2};
(* Represent dM/dt using dq as generalized velocities *)
(* Here one would substitute specific M and C, then evaluate x.S.x // Simplify *)

      

For a concrete small robot (e.g. a planar 2-DOF arm), you can define explicit expressions for \( \mathbf{M}(\mathbf{q}) \), \( \mathbf{C}(\mathbf{q},\dot{\mathbf{q}}) \), \( \mathbf{g}(\mathbf{q}) \), and \( U(\mathbf{q}) \), then symbolically verify that \( \mathbf{x}^\top(\dot{\mathbf{M}} - 2\mathbf{C})\mathbf{x} = 0 \) and \( \dot{H} = \dot{\mathbf{q}}^\top \boldsymbol{\tau} \).

In this lesson we connected the structure of the Lagrange–Euler equations with fundamental energetic properties of manipulator dynamics. Starting from the kinetic energy representation and the gravity potential, we showed that:

• The inertia matrix \( \mathbf{M}(\mathbf{q}) \) is symmetric positive definite.
• The skew-symmetry of \( \dot{\mathbf{M}}(\mathbf{q}) - 2\mathbf{C}(\mathbf{q},\dot{\mathbf{q}}) \) implies exact power balance \( \dot{H} = \dot{\mathbf{q}}^\top \boldsymbol{\tau} \).
• Mechanical energy provides a natural storage function certifying passivity from joint torques to joint velocities, with viscous damping yielding output strict passivity.
• Practical energy consistency and passivity checks can be implemented numerically across different programming environments by checking matrix properties, skew-symmetry conditions, and simulated energy trajectories.

These properties are crucial for trusting derived dynamic models and will be used in the next lesson, where we apply them to concrete 2-DOF and 3-DOF manipulator examples.