In classical rigid-body manipulation, a robot interacts with objects whose configuration can be described by a small number of generalized coordinates. For example, with \( m \) rigid objects and a robot with \( n \) joints, a typical configuration space is

\[ \mathcal{Q}_{\text{rigid}} \cong SE(3)^m \times \mathbb{R}^n, \]

which is a smooth manifold of finite dimension \( \dim(\mathcal{Q}_{\text{rigid}}) = 6m + n \). Algorithms from earlier chapters (graph search, sampling-based planning, trajectory optimization) crucially rely on this finite-dimensional structure.

By contrast, a deformable object (cloth, rope, soft bag) has a configuration that is a field over a continuous domain. Intuitively, there are infinitely many degrees of freedom: each point of the material can move in space. Even when we discretize, the resulting configuration can easily be hundreds or thousands of dimensions, with strong coupling and numerical stiffness.

Deformable manipulation is harder than rigid manipulation for at least four tightly coupled reasons:

• High-dimensional state: the configuration space is infinite-dimensional in principle, and very high-dimensional after discretization.
• Complex dynamics and stiffness: elastic forces lead to stiff differential equations that are hard to integrate stably.
• Contact and underactuation: the robot touches a deformable body only at a few contact points, yet must control its global shape.
• Partial observability: sensors see only parts of the object (occlusion, self-contact), making state estimation ill-posed.

flowchart LR
  R["Rigid manipulation \n(few DOFs)"] --> D["Deformable manipulation \n(many DOFs)"]
  D --> S1["High-dimensional state"]
  D --> S2["Complex contact constraints"]
  D --> S3["Partial and noisy sensing"]
  D --> S4["Numerical stiffness in simulation"]
        

Consider a deformable body occupying a reference domain \( \Omega \subset \mathbb{R}^3 \) in material (Lagrangian) coordinates. A configuration at time \( t \) is a mapping

\[ x(\cdot, t) : \Omega \to \mathbb{R}^3, \quad X \mapsto x(X,t), \]

where \( X \in \Omega \) is a material point and \( x(X,t) \) is its position in space. The configuration space of the deformable body is then the function space

\[ \mathcal{Q}_{\text{def}} = \{\, x : \Omega \to \mathbb{R}^3 \mid x \text{ satisfies boundary conditions and regularity} \,\}. \]

This is infinite-dimensional: specifying a configuration means specifying a vector in \( \mathbb{R}^3 \) for each material point.

In practice, we approximate \( \mathcal{Q}_{\text{def}} \) by a finite-dimensional subspace using discretization (particles, mass–spring systems, finite elements). For example, a simple rope model samples \( N \) material points along the arc-length coordinate with positions \( x_i \in \mathbb{R}^3 \) for \( i=1,\dots,N \). The stacked configuration vector is

\[ q = \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_N \end{bmatrix} \in \mathbb{R}^{3N}, \]

so the discrete configuration space is \( \mathbb{R}^{3N} \). As \( N \) grows, this finite-dimensional approximation becomes more accurate but also more challenging for planning and control algorithms.

A common energy for a mass–spring rope is the elastic potential

\[ U(q) = \frac{1}{2} \sum_{i=1}^{N-1} k_i \Big( \|x_{i+1} - x_i\| - \ell_i \Big)^2, \]

where \( k_i > 0 \) is the stiffness of the spring between nodes \( i \) and \( i+1 \), and \( \ell_i \) is its rest length. Even this simple quadratic energy introduces long-range coupling: moving one node affects forces at neighboring nodes through the gradient \( \nabla U(q) \).

For a rigid manipulator, you have already seen the familiar dynamics

\[ M_r(q_r)\ddot{q}_r + C_r(q_r,\dot{q}_r)\dot{q}_r + g_r(q_r) = \boldsymbol{\tau}, \]

where \( q_r \in \mathbb{R}^n \) and \( M_r \) is positive definite. For a discrete deformable body with configuration \( q \in \mathbb{R}^{3N} \), a typical semi-discrete model is

\[ M \ddot{q} + D \dot{q} + \nabla U(q) = f_{\text{ext}}(q,\dot{q},t), \]

where:

• \( M \in \mathbb{R}^{3N \times 3N} \) is the mass matrix (often diagonal for particles).
• \( D \) is a damping matrix modeling internal friction.
• \( U(q) \) is the elastic potential (e.g., springs, FEM).
• \( f_{\text{ext}} \) contains gravity and contact forces from the robot and environment.

Linearizing the forces around an equilibrium \( q^\star \) gives

\[ M \ddot{\delta q} + D \dot{\delta q} + K \,\delta q \approx 0, \quad K = \nabla^2 U(q^\star), \]

where \( K \) is the stiffness matrix. The eigenvalues of \( M^{-1}K \) determine oscillation frequencies. If the largest eigenvalue \( \lambda_{\max} \) is large, the system is stiff: some modes oscillate very quickly and force tiny time steps in explicit integration schemes.

For a simple explicit (forward) Euler time integration of \( M \ddot{q} + K q = 0 \), a standard stability bound is

\[ \Delta t \le \frac{2}{\sqrt{\lambda_{\max}(M^{-1}K)}}. \]

As we refine the discretization (increase \( N \) or stiffness), \( \lambda_{\max} \) tends to grow, shrinking the stable step size \( \Delta t \). This creates a fundamental trade-off:

• Coarse models: easier to simulate and plan with, but inaccurate shapes.
• Fine models: accurate, but numerically stiff and costly.

This trade-off is one of the core reasons deformable planning often uses reduced-order models or very coarse approximations in practice.

A rigid gripper typically makes contact with a deformable object at a small number of points or patches. Let \( \theta \in \mathbb{R}^m \) denote the robot joint coordinates and \( q \in \mathbb{R}^{3N} \) the deformable configuration. Locally, we may write the coupling as

\[ q = \Phi(\theta), \quad \dot{q} = J(\theta)\dot{\theta}, \quad J(\theta) = \frac{\partial \Phi}{\partial \theta}(\theta). \]

Since the robot controls only \( m \) joints, \( \operatorname{rank}(J(\theta)) \le m \). However, \( q \) has \( 3N \) components, with \( 3N \gg m \) in realistic models. This means:

• Only a low-dimensional subspace of deformations is directly actuated.
• Many internal deformation modes lie in the nullspace of \( J(\theta) \), and cannot be directly controlled.

At contact, we must also obey unilateral (non-penetration) constraints and friction laws. For a normal gap function \( g(q) \) between the deformable surface and an environmental obstacle, we have complementarity conditions

\[ g(q) \ge 0, \quad \lambda_n \ge 0, \quad g(q)\,\lambda_n = 0, \]

where \( \lambda_n \) is the normal contact force multiplier. These relations define a feasible set for \( (q,\lambda_n) \) that is highly nonconvex. When combined with friction cones and the high-dimensional state, this makes planning in the joint space \( \theta \) extremely challenging: small changes in motion can cause large, discontinuous changes in contact configuration and object shape.

In rigid-body manipulation, we often estimate a small pose vector (e.g., a single transform in \( SE(3) \)) from images or depth data. For deformables, we must estimate the entire high-dimensional state \( q \) from limited measurements \( y \). A generic observation model is

\[ y = h(q) + \eta, \]

where \( \eta \) is sensor noise. Linearizing around a reference \( q^\star \) gives

\[ y \approx h(q^\star) + H(q^\star)(q - q^\star), \quad H(q^\star) = \frac{\partial h}{\partial q}\bigg|_{q^\star}. \]

If we collect \( m_{\text{meas}} \) scalar measurements and \( q \in \mathbb{R}^{3N} \), the Jacobian \( H(q^\star) \in \mathbb{R}^{m_{\text{meas}} \times 3N} \) typically has \( m_{\text{meas}} \ll 3N \). Then the linear system for \( q \) is severely underdetermined, and many different internal shapes explain the same sensor readings.

In practice, this forces us to:

• Introduce priors on plausible shapes (smoothness, low rank, low curvature).
• Use low-dimensional embeddings (e.g., modal coordinates) for \( q \).
• Accept that the robot may never know the true full state, only a coarse approximation relevant for the task.

The combination of underactuation (Section 4) and partial observability means that many classical control design tools relying on full-state feedback and exact models become fragile or inapplicable.

To concretize the above ideas, we simulate a planar rope as a mass–spring chain in Python. One end is fixed, and the other end is directly driven by the robot gripper along a sinusoidal trajectory. The rope is modeled by \( N \) point masses coupled by springs and linear damping. Note how small the time step must be to keep the simulation stable.

import numpy as np

# Number of nodes in the rope
N = 10
mass = 0.05
k_spring = 50.0
damping = 0.1
g = 9.81

# Small time step due to stiffness
dt = 5e-4
steps = 20000

# Positions and velocities in 2D
x = np.zeros((N, 2))
v = np.zeros_like(x)

# Initialize rope horizontally
rest_length = 0.05
for i in range(N):
    x[i, 0] = rest_length * i
    x[i, 1] = 0.0

rest_lengths = np.full(N - 1, rest_length)

def compute_forces(x, v):
    f = np.zeros_like(x)

    # Internal springs between neighbors
    for i in range(N - 1):
        dx = x[i + 1] - x[i]
        L = np.linalg.norm(dx)
        if L > 1e-6:
            direction = dx / L
            f_spring = k_spring * (L - rest_lengths[i]) * direction
            f[i] += f_spring
            f[i + 1] -= f_spring

    # Damping and gravity
    f -= damping * v
    f[:, 1] -= g * mass
    return f

for step in range(steps):
    # External forces
    f = compute_forces(x, v)

    # Fix the first node (anchored to the table)
    f[0, :] = 0.0
    v[0, :] = 0.0
    x[0, :] = np.array([0.0, 0.0])

    # Drive the last node along a simple trajectory (robot gripper)
    t = step * dt
    target = np.array([rest_length * (N - 1),
                       0.1 * np.sin(2.0 * np.pi * t)])
    x[-1, :] = target
    v[-1, :] = np.array([0.0, 0.0])

    # Semi-implicit (symplectic) Euler integration
    a = f / mass
    v += dt * a
    x += dt * v

# In practice, this core integration loop would be embedded in a robotics
# simulator (e.g., PyBullet, MuJoCo), and coupled to robot joint-space control.
      

If we increase \( k_{\text{spring}} \) or \( N \) without adjusting \( dt \), the simulation quickly becomes unstable, numerically illustrating the stiffness issues discussed in Section 3.

In C++, the Eigen library is a natural choice for vector and matrix operations. The following minimal example mirrors the Python code.

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

int main() {
    const int N = 10;
    const double mass = 0.05;
    const double k_spring = 50.0;
    const double damping = 0.1;
    const double g = 9.81;
    const double dt = 5e-4;
    const int steps = 20000;

    std::vector<Eigen::Vector2d> x(N), v(N);
    for (int i = 0; i < N; ++i) {
        x[i] = Eigen::Vector2d(0.05 * i, 0.0);
        v[i].setZero();
    }
    std::vector<double> rest_lengths(N - 1, 0.05);

    auto compute_forces =
        [&](const std::vector<Eigen::Vector2d> &x,
              const std::vector<Eigen::Vector2d> &v) {
            std::vector<Eigen::Vector2d> f(N, Eigen::Vector2d::Zero());
            for (int i = 0; i < N - 1; ++i) {
                Eigen::Vector2d dx = x[i + 1] - x[i];
                double L = dx.norm();
                if (L > 1e-6) {
                    Eigen::Vector2d dir = dx / L;
                    Eigen::Vector2d fs = k_spring * (L - rest_lengths[i]) * dir;
                    f[i] += fs;
                    f[i + 1] -= fs;
                }
            }
            for (int i = 0; i < N; ++i) {
                f[i] -= damping * v[i];
                f[i].y() -= g * mass;
            }
            return f;
        };

    for (int step = 0; step < steps; ++step) {
        auto f = compute_forces(x, v);

        // Anchor first node
        f[0].setZero();
        v[0].setZero();
        x[0] = Eigen::Vector2d(0.0, 0.0);

        // Drive last node (gripper)
        double t = step * dt;
        Eigen::Vector2d target(0.05 * (N - 1),
                               0.1 * std::sin(2.0 * M_PI * t));
        x[N - 1] = target;
        v[N - 1].setZero();

        for (int i = 0; i < N; ++i) {
            Eigen::Vector2d a = f[i] / mass;
            v[i] += dt * a;
            x[i] += dt * v[i];
        }
    }

    std::cout << "Simulation finished." << std::endl;
    return 0;
}
      

A simple Java version can be implemented with primitive arrays or a linear algebra library (e.g., EJML). Below is a bare-bones implementation using double[][].

public class RopeSimulator {
    private final int N = 10;
    private final double mass = 0.05;
    private final double kSpring = 50.0;
    private final double damping = 0.1;
    private final double g = 9.81;
    private final double dt = 5e-4;

    private final double[][] x = new double[N][2];
    private final double[][] v = new double[N][2];
    private final double[] restLengths = new double[N - 1];

    public RopeSimulator() {
        double rest = 0.05;
        for (int i = 0; i < N; ++i) {
            x[i][0] = rest * i;
            x[i][1] = 0.0;
            v[i][0] = 0.0;
            v[i][1] = 0.0;
        }
        for (int i = 0; i < N - 1; ++i) {
            restLengths[i] = rest;
        }
    }

    private void computeForces(double[][] f) {
        for (int i = 0; i < N; ++i) {
            f[i][0] = 0.0;
            f[i][1] = 0.0;
        }

        // Springs
        for (int i = 0; i < N - 1; ++i) {
            double dx = x[i + 1][0] - x[i][0];
            double dy = x[i + 1][1] - x[i][1];
            double L = Math.sqrt(dx * dx + dy * dy);
            if (L > 1e-6) {
                double dirx = dx / L;
                double diry = dy / L;
                double fs = kSpring * (L - restLengths[i]);
                double fx = fs * dirx;
                double fy = fs * diry;
                f[i][0] += fx;
                f[i][1] += fy;
                f[i + 1][0] -= fx;
                f[i + 1][1] -= fy;
            }
        }

        // Damping and gravity
        for (int i = 0; i < N; ++i) {
            f[i][0] -= damping * v[i][0];
            f[i][1] -= damping * v[i][1];
            f[i][1] -= g * mass;
        }
    }

    public void step(int stepIndex) {
        double[][] f = new double[N][2];
        computeForces(f);

        // Anchor first node
        f[0][0] = 0.0;
        f[0][1] = 0.0;
        v[0][0] = 0.0;
        v[0][1] = 0.0;
        x[0][0] = 0.0;
        x[0][1] = 0.0;

        // Drive last node
        double t = stepIndex * dt;
        x[N - 1][0] = 0.05 * (N - 1);
        x[N - 1][1] = 0.1 * Math.sin(2.0 * Math.PI * t);
        v[N - 1][0] = 0.0;
        v[N - 1][1] = 0.0;

        for (int i = 0; i < N; ++i) {
            double ax = f[i][0] / mass;
            double ay = f[i][1] / mass;
            v[i][0] += dt * ax;
            v[i][1] += dt * ay;
            x[i][0] += dt * v[i][0];
            x[i][1] += dt * v[i][1];
        }
    }

    public static void main(String[] args) {
        RopeSimulator sim = new RopeSimulator();
        int steps = 20000;
        for (int step = 0; step < steps; ++step) {
            sim.step(step);
        }
        System.out.println("Simulation finished.");
    }
}
      

In a robotics system, the positions of the end node would be linked to the forward kinematics of the robot gripper, and the entire simulation would be embedded in a control loop or planning algorithm.

MATLAB and Simulink are well-suited for prototyping deformable models as state-space systems. Below is a MATLAB script that builds the mass and stiffness matrices for a 1D chain and simulates it with a simple explicit integrator. The same \( M \), \( K \), and damping can be implemented in Simulink using a State-Space block.

N = 10;
m = 0.05;
k_spring = 50;
c_damp = 0.1;
g = 9.81;

% Mass and stiffness for 1D vertical motion per node (y-coordinate only)
M = m * eye(N);
K = zeros(N);
for i = 1:N-1
    e = zeros(N,1);
    e(i) = 1; e(i+1) = -1;
    K = K + k_spring * (e * e');
end

% Damping proportional to mass
C = c_damp * eye(N);

dt = 5e-4;
steps = 20000;
y = zeros(N,1);
v = zeros(N,1);

% Initialization (rope horizontal, small sag)
for i = 1:N
    y(i) = 0.0;
end

for step = 1:steps
    % Gravity
    f_ext = -g * m * ones(N,1);

    % Anchor first node
    f_ext(1) = 0;
    y(1) = 0;
    v(1) = 0;

    % Drive last node (robot gripper)
    t = step * dt;
    y(N) = 0.1 * sin(2*pi*t);
    v(N) = 0;

    % Internal elastic and damping forces: M * a = f_ext - C*v - K*y
    rhs = f_ext - C * v - K * y;
    a = M \ rhs;

    % Semi-implicit Euler
    v = v + dt * a;
    y = y + dt * v;
end

% In Simulink, put M, C, K into a second-order state-space block and
% connect the last node to the robot gripper trajectory as an input.
      

In Wolfram Mathematica, we can use NDSolve to integrate the ODEs of a small chain. For illustration, consider a 2-mass system:

m = 0.05;
k = 50;
c = 0.1;
g = 9.81;

(* y1[t], y2[t] are vertical positions of two masses;
   y1[t] is fixed at zero to model an anchor. *)
eqs = {
  y1[t] == 0,
  m y2''[t] + c y2'[t] + k (y2[t] - y1[t]) == -m g
};

ics = {
  y2[0] == 0,
  y2'[0] == 0
};

sol = NDSolve[
  Join[eqs, ics],
  {y1, y2},
  {t, 0, 2}
];

Plot[Evaluate[y2[t] /. sol], {t, 0, 2},
  PlotLabel -> "Second mass trajectory"]
      

Extending this construction to longer chains yields mass–spring ropes analogous to the Python/C++/Java examples. One can then build control laws at the level of generalized coordinates \( q \).

The following diagram summarizes the typical workflow for using deformable models in robotic manipulation: start from continuum mechanics, discretize to obtain a tractable state representation, then design planning and control algorithms on this reduced model and close the loop with sensing.

flowchart TD
  C["Continuum model (PDEs on domain)"] --> D["Spatial discretization (particles, FEM)"]
  D --> T["Time integration and simulation"]
  T --> P["Planning and control on discrete state"]
  P --> E["Execution on physical robot"]
  E --> F["Feedback and model update"]
        

Each arrow hides nontrivial numerical analysis and approximation choices, and errors introduced in early stages (e.g., overly coarse meshes) can fundamentally limit what manipulation tasks are achievable.

A simple linearized model around a nominal configuration \( q^\star \) can be written as

\[ \delta \dot{q} = A \,\delta q + B \,\delta \theta, \]

where \( \delta q \in \mathbb{R}^{3N} \) and \( \delta \theta \in \mathbb{R}^m \). The input matrix \( B \) is related to the Jacobian \( J(\theta^\star) \) in Section 4. Because \( m \ll 3N \), the rank of \( B \) is at most \( m \). The controllable subspace of the pair \( (A,B) \) therefore has dimension at most \( m \) (and often less, depending on \( A \)).

This implies that many internal deformation modes (eigenvectors of \( A \) or \( M^{-1}K \)) are not directly controllable. In practice, robots exploit:

• Passive dynamics: allowing the rope or cloth to settle into low-energy shapes.
• Environmental constraints: dragging a cloth against a table to fold it.
• Heuristic strategies: aiming for task-level goals (e.g., “endpoints aligned”) rather than full-state tracking.

Even when full controllability fails, task-level controllability can hold: for instance, the robot might not control every point of a cable, but can reliably place its endpoints in desired locations.

Problem 1 (Dimension of a Discrete Rope Configuration Space): A rope in 3D is modeled as \( N \) point masses with positions \( x_i \in \mathbb{R}^3 \), \( i=1,\dots,N \). One end (\( x_1 \)) is fixed at the origin. What is the dimension of the configuration space of the remaining masses? How does this compare with a single rigid body in \( SE(3) \)?

Solution: Without constraints, the stacked vector \( q = [x_1^\top,\dots,x_N^\top]^\top \) lies in \( \mathbb{R}^{3N} \). Fixing \( x_1 = 0 \) removes three degrees of freedom, so the configuration space is isomorphic to \( \mathbb{R}^{3N-3} \), of dimension \( 3N-3 \). For a single rigid body in 3D, the configuration space is \( SE(3) \) with dimension 6. Thus, even for modest \( N \), the rope’s configuration space is far higher dimensional than that of a rigid body.

Problem 2 (Stability Constraint on Time Step): Suppose a 1D mass–spring chain has mass matrix \( M = m I \) and stiffness matrix \( K \) whose largest eigenvalue satisfies \( \lambda_{\max}(M^{-1}K) = \omega_{\max}^2 \). Show that the explicit Euler scheme \( q_{k+1} = q_k + \Delta t \,\dot{q}_k \), \( \dot{q}_{k+1} = \dot{q}_k - \Delta t \,M^{-1}K q_k \) is stable only if \( \Delta t \le 2 / \omega_{\max} \).

Solution: For a single mode with amplitude \( z \), the linearized system is

\[ \ddot{z} + \omega_{\max}^2 z = 0. \]

Introducing \( v = \dot{z} \), explicit Euler yields the update

\[ \begin{bmatrix} z_{k+1} \\ v_{k+1} \end{bmatrix} = \begin{bmatrix} 1 & \Delta t \\ -\omega_{\max}^2 \Delta t & 1 \end{bmatrix} \begin{bmatrix} z_k \\ v_k \end{bmatrix}. \]

The eigenvalues \( \rho \) of the update matrix satisfy

\[ \rho^2 - 2\rho + (1 + \omega_{\max}^2 \Delta t^2) = 0. \]

Stability of the discrete-time system requires \( |\rho| \le 1 \). Solving the quadratic shows \( |\rho| > 1 \) whenever \( \omega_{\max}\Delta t > 2 \). Therefore the necessary stability condition is \( \Delta t \le 2/\omega_{\max} \), consistent with the bound in Section 3.

Problem 3 (Underactuated Nullspace): Let \( q \in \mathbb{R}^{3N} \) and \( \theta \in \mathbb{R}^m \). Suppose locally \( q = \Phi(\theta) \) and \( J(\theta) = \partial \Phi / \partial \theta \) has rank \( m \). Show that the dimension of the nullspace of \( J(\theta)^\top \) is at least \( 3N - m \), and interpret this physically.

Solution: Since \( J(\theta) \in \mathbb{R}^{3N \times m} \) has rank \( m \), the rank–nullity theorem implies that \( \dim(\ker J(\theta)) = 3N - m \). The matrix \( J(\theta)^\top \) has the same rank \( m \), so

\[ \dim(\ker J(\theta)^\top) = 3N - m. \]

Vectors in \( \ker J(\theta)^\top \) correspond to generalized forces on the deformable object that do not map to any joint torques. Physically, they represent deformation modes that cannot be directly generated (or sensed) by the robot’s actuators and sensors, highlighting strong underactuation.

Problem 4 (Partial Observability): Assume the linearized observation model \( y = H q + \eta \), with \( H \in \mathbb{R}^{m_{\text{meas}} \times 3N} \) and \( m_{\text{meas}} < 3N \). Show that the least-squares estimate \( \hat{q}_{\text{LS}} \) is not unique unless additional priors are added.

Solution: The least-squares objective is

\[ \min_q \| H q - y \|^2. \]

The normal equations are \( H^\top H q = H^\top y \). If \( m_{\text{meas}} < 3N \), then \( \operatorname{rank}(H) \le m_{\text{meas}} < 3N \), so \( H^\top H \) is singular and has a nontrivial nullspace. Therefore, any \( q_0 \) satisfying the normal equations generates a whole affine subspace \( \{ q_0 + z \mid z \in \ker(H) \} \) of minimizers. Regularization (e.g., penalizing \( \|q\|^2 \) or curvature) or low-dimensional parameterizations are required to obtain a unique, physically plausible estimate.

In this lesson we analyzed why deformable object manipulation is substantially harder than rigid manipulation for a rigid robot. The key difficulties are: (i) infinite-dimensional configuration spaces that must be approximated by high-dimensional discrete models; (ii) stiff dynamics that constrain time stepping and demand careful numerical methods; (iii) underactuation and nonconvex contact constraints that limit controllability and complicate planning; and (iv) partial observability that makes state estimation ill-posed without strong priors. We grounded these ideas in a simple mass–spring rope model implemented in Python, C++, Java, MATLAB/Simulink, and Mathematica, setting the stage for subsequent lessons on state representations, planning under uncertainty, and learning-based methods for deformable manipulation.