A robot interacting with its environment is naturally modeled as a hybrid system with continuous-time dynamics inside each mode and discrete jumps when contacts appear or disappear. Typical modes for a manipulator with a single unilateral contact are:

• Free motion: no active contact forces; only gravity, Coriolis, and actuation act on the robot.
• Continuous contact: normal distance constraint is active, with continuous contact forces enforcing non-penetration.
• Impact event: an instantaneous velocity jump when the end-effector collides with the environment.

Let \( q \in \mathbb{R}^n \) be joint coordinates and \( v = \dot{q} \) joint velocities. In all modes, the robot obeys the standard manipulator dynamics

\[ \mathbf{M}(q)\dot{v} + \mathbf{h}(q,v) = \mathbf{S} \tau + \mathbf{J}_c(q)^\top f_c , \]

where \( \mathbf{M}(q) \) is the inertia matrix, \( \mathbf{h}(q,v) \) collects Coriolis and gravity, \( \mathbf{S} \) selects actuated joints, \( \tau \) are joint torques, and \( \mathbf{J}_c(q) \) is the Jacobian of the contact point(s). The contact force \( f_c \) is mode dependent: in free motion it is zero; in continuous contact it is determined by enforcing kinematic constraints; during impact it is replaced by a contact impulse.

We encapsulate the main modes and their transitions in the following finite-state representation.

flowchart TD
  F["Free motion"] -->|"approach surface (gap = 0 and normal vel <= 0)"| I["Impact event"]
  I -->|"velocity reset"| C["Continuous contact"]
  C -->|"separation (normal force = 0 and normal vel >= 0)"| F
  C -->|"slip / stick changes"| C
        

The controller must (i) respect these modes and transitions by enforcing non-penetration and (ii) shape motion so that impacts are either benign (low-energy) or avoided when necessary.

Consider a manipulator whose end-effector contacts a rigid surface with a normal direction given by a contact Jacobian row \( \mathbf{J}_n(q) \in \mathbb{R}^{1 \times n} \). At the instant of contact, the configuration \( q \) is assumed continuous, but the velocity \( v \) may undergo a discontinuous jump:

\[ v^+ = v^- + \Delta v, \]

where \( v^- \) is the pre-impact velocity, \( v^+ \) is the post-impact velocity, and \( \Delta v \) is determined by the impulse-momentum balance.

2.1 Impulse-Momentum Balance

Let \( \Lambda_n \) be the normal impulse (integral of normal contact force over the infinitesimal impact time). For simplicity we ignore tangential frictional impulses here. The joint-space impulse is \( \mathbf{S}_J = \mathbf{J}_n(q)^\top \Lambda_n \). Over the short impact interval \( [t^-, t^+] \), we neglect actuation and smooth forces, so the only significant effect is the impulse:

\[ \mathbf{M}(q)\left( v^+ - v^- \right) = \mathbf{J}_n(q)^\top \Lambda_n . \]

We enforce that after impact, the end-effector normal velocity satisfies a restitution law. Let the pre- and post-impact normal velocities be \( \dot{x}_n^- = \mathbf{J}_n(q)v^- \) and \( \dot{x}_n^+ = \mathbf{J}_n(q)v^+ \). A standard (one-dimensional) Newton restitution model is

\[ \dot{x}_n^+ = -e\, \dot{x}_n^- , \]

where \( 0 \le e \le 1 \) is the coefficient of restitution (\( e = 0 \) perfectly inelastic, \( e = 1 \) perfectly elastic).

2.2 Solving for Post-Impact Velocity

Substitute the impulse-momentum equation into the restitution law. From \( v^+ = v^- + \mathbf{M}(q)^{-1}\mathbf{J}_n(q)^\top \Lambda_n \), we get

\[ \dot{x}_n^+ = \mathbf{J}_n(q)v^+ = \mathbf{J}_n(q)v^- + \mathbf{J}_n(q)\mathbf{M}(q)^{-1}\mathbf{J}_n(q)^\top \Lambda_n . \]

Imposing \( \dot{x}_n^+ = -e\,\dot{x}_n^- \) yields

\[ -e\,\dot{x}_n^- = \dot{x}_n^- + \mathbf{J}_n(q)\mathbf{M}(q)^{-1} \mathbf{J}_n(q)^\top \Lambda_n . \]

Solve for \( \Lambda_n \):

\[ \Lambda_n = -\frac{(1+e)\,\dot{x}_n^-} {\mathbf{J}_n(q)\mathbf{M}(q)^{-1}\mathbf{J}_n(q)^\top}. \]

Provided the pre-impact normal velocity satisfies \( \dot{x}_n^- < 0 \) (approaching the surface), the numerator is positive and we obtain a positive normal impulse, consistent with repulsive contact.

Substituting back gives the post-impact velocity:

\[ v^+ = v^- - \mathbf{M}(q)^{-1}\mathbf{J}_n(q)^\top \frac{(1+e)\,\dot{x}_n^-} {\mathbf{J}_n(q)\mathbf{M}(q)^{-1}\mathbf{J}_n(q)^\top}. \]

2.3 Energy Loss at Impact (Sketch)

The kinetic energy before impact is \( K^- = \tfrac{1}{2} {v^-}^\top \mathbf{M}(q)v^- \), and after impact \( K^+ = \tfrac{1}{2} {v^+}^\top \mathbf{M}(q)v^+ \). One can show (by substituting the expression for \( v^+ \) and simplifying) that

\[ K^+ - K^- = -\frac{1-e^2}{2} \frac{\left(\dot{x}_n^-\right)^2} {\mathbf{J}_n(q)\mathbf{M}(q)^{-1}\mathbf{J}_n(q)^\top} \le 0, \]

that is, kinetic energy never increases at impact when \( 0 \le e \le 1 \). This fact is important for later Lyapunov-based stability analyses of hybrid robot controllers.

We formalize contact transitions using a hybrid automaton. The state is the pair \( x = (q,v) \) and the discrete mode \( \sigma \in \{\text{free}, \text{contact}\} \). For each mode, we have ordinary differential equations (ODEs) describing continuous evolution and guards specifying when jumps occur.

\[ \dot{x} = f_\sigma(x,\tau) = \begin{bmatrix} v \\ \mathbf{M}(q)^{-1} \big(\mathbf{S} \tau - \mathbf{h}(q,v) + \mathbf{J}_c(q)^\top f_c(x,\sigma)\big) \end{bmatrix}. \]

• In free mode, \( f_c = 0 \) and the state evolves without contact forces until the guard \( g_\text{impact}(x) = \phi(q) = 0 \) and \( \dot{\phi}(q,v) < 0 \) is satisfied (here \( \phi(q) \) is the signed distance to the surface).
• At the guard, a jump is triggered: \( x^+ = \Delta(x^-) \), where \( \Delta \) applies the impact law from the previous section (velocity reset, configuration unchanged).
• In contact mode, we enforce the holonomic constraint \( \phi(q)=0 \) and \( \dot{\phi}(q,v)=0 \), which determines \( f_c \). When the normal force reaches zero and the normal acceleration becomes positive, the contact can lift off.

Mathematically, this leads to complementarity conditions such as

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

where \( \lambda_n \) is the normal component of the contact force in continuous mode. Solving these relations together with the dynamics provides consistent contact-aware motion.

Many joint- or task-space controllers are derived under the assumption of smooth dynamics, but contact transitions violate this assumption. We outline simple impact-aware modifications for a model-based controller of the form

\[ \tau = \mathbf{M}(q)\nu + \mathbf{h}(q,v) - \mathbf{J}_c(q)^\top f_c^\text{des}, \]

where \( \nu \) is a desired joint acceleration obtained from tracking errors and \( f_c^\text{des} \) is a desired contact force (possibly zero in free motion).

4.1 Pre-Impact Velocity Shaping

Assume that we know the nominal time or configuration where the end-effector will hit a surface. We can design a pre-contact phase where the task-space normal velocity is regulated to a desired impact value \( \dot{x}_n^\text{des} \):

\[ \dot{x}_n = \mathbf{J}_n(q)v,\quad \ddot{x}_n = \mathbf{J}_n(q)\dot{v} + \dot{\mathbf{J}}_n(q,v)v. \]

Using inverse dynamics and a simple feedback term,

\[ \ddot{x}_n^\text{cmd} = -k_n\big(\dot{x}_n - \dot{x}_n^\text{des}\big), \]

and embedding this into a task-space control law, we can ensure that the impact speed is small, thereby reducing impact energy. Practically, this is combined with tangential tracking objectives in a weighted fashion.

4.2 Post-Impact Mode Switching and Damping Injection

Immediately after impact, the controller should:

1. Switch from a free-motion control law to a constraint-consistent control law that enforces \( \phi(q)=0 \).
2. Increase damping gains along the normal direction to dissipate residual kinetic energy.

In a joint-space controller, this can be achieved by defining a contact-aware gain matrix \( \mathbf{K}_d^\sigma \) that depends on the mode \( \sigma \):

\[ \tau = \mathbf{M}(q)\ddot{q}^\text{des} + \mathbf{h}(q,v) - \mathbf{K}_p (q - q^\text{des}) - \mathbf{K}_d^\sigma (v - v^\text{des}). \]

For example, in free mode one may choose uniform diagonal damping, whereas in contact mode one projects a larger damping onto the normal subspace spanned by \( \mathbf{J}_n(q) \), leaving tangential directions less affected to maintain tracking performance.

A high-level algorithmic view of an impact-aware control loop is shown below.

flowchart TD
  S["Read q, v at time t"] --> U["Compute control torque tau(q, v, sigma)"]
  U --> INT["Integrate dynamics until next step or event"]
  INT --> G{"Guard fired \n(impact or liftoff)?"}
  G -->|no| S
  G -->|yes| R["Apply reset map (impact law or mode change)"]
  R --> M["Update mode sigma and contact model"]
  M --> S
        

We implement a simple 1-DOF robot (mass on a vertical line) with a rigid stop at \( q = 0 \). The state is \( x = (q,v) \), with gravity \( g > 0 \) pulling downward (\( q \) decreases). An impact occurs when \( q = 0 \) and \( v < 0 \).

This prototype demonstrates:

• Contact event detection.
• Velocity reset using a coefficient of restitution.
• Simple PD gravity-compensated control above the surface.

import numpy as np

class OneDOFImpactSystem:
    def __init__(self, m=1.0, g=9.81, e=0.0):
        self.m = m      # mass
        self.g = g      # gravity
        self.e = e      # restitution (0 inelastic, 1 elastic)

    def dynamics(self, t, x, tau):
        """
        Continuous dynamics between impacts.
        x = [q, v]. tau is control torque (force).
        """
        q, v = x
        # mass * acceleration = tau - m*g
        a = (tau - self.m * self.g) / self.m
        return np.array([v, a])

    def impact_event(self, t, x):
        """
        Guard function: impact when q = 0 and v < 0.
        We only encode the surface condition here; direction handled separately.
        """
        q, _ = x
        return q  # zero when contact with the stop

    def handle_impact(self, x):
        """
        Apply Newton impact law with coefficient of restitution e.
        """
        q, v = x
        if q == 0.0 and v < 0.0:
            v_plus = -self.e * v  # normal velocity reversed and scaled
            return np.array([q, v_plus])
        else:
            return x

def pd_controller(q, v, q_des, v_des, kp=50.0, kd=10.0, g=9.81, m=1.0):
    """
    Simple PD with gravity compensation for a vertical 1-DOF system.
    """
    # gravity compensation
    tau_g = m * g
    # PD term
    tau_pd = -kp * (q - q_des) - kd * (v - v_des)
    return tau_g + tau_pd

def simulate(system, x0, t0, tf, dt, q_des):
    t = t0
    x = np.array(x0, dtype=float)
    ts = [t]
    xs = [x.copy()]

    while t < tf:
        q, v = x
        # Mode: free or resting on the stop
        if q > 0.0 or v > 0.0:
            # free motion (above or separating from the stop)
            tau = pd_controller(q, v, q_des, 0.0, g=system.g, m=system.m)
        else:
            # we are on the stop and not separating (q == 0 and v == 0)
            # hold position with PD to avoid drift
            tau = 0.0

        # integrate one Euler step (for simplicity)
        dx = system.dynamics(t, x, tau)
        x_new = x + dt * dx

        # detect surface crossing q = 0
        if x[0] > 0.0 and x_new[0] <= 0.0:
            # approximate impact at intermediate time
            # move to surface
            alpha = x[0] / (x[0] - x_new[0] + 1e-12)
            x_impact = x + alpha * (x_new - x)

            # apply impact law
            x_after = system.handle_impact(x_impact)

            # continue with remaining time in the step if needed
            remaining = (1.0 - alpha) * dt
            dx_after = system.dynamics(t, x_after, tau)
            x_new = x_after + remaining * dx_after

        t += dt
        x = x_new
        ts.append(t)
        xs.append(x.copy())

    return np.array(ts), np.stack(xs, axis=0)

if __name__ == "__main__":
    sys = OneDOFImpactSystem(m=1.0, g=9.81, e=0.2)
    x0 = [0.5, -1.0]  # start above the stop, moving downward
    t, x = simulate(sys, x0, t0=0.0, tf=3.0, dt=1e-3, q_des=0.2)
    # Here you can plot q(t) and v(t) to inspect impact behaviour.
      

In more complex manipulators one would use robotics libraries such as pinocchio, pybullet, or ROS-based ros_control to obtain \( \mathbf{M}(q) \), \( \mathbf{h}(q,v) \), and contact Jacobians, then embed similar impact handling in their simulation or control loop.

6.1 C++ Skeleton (Using a Generic Rigid-Body Library)

The following C++ snippet outlines how one might structure a control loop with impact-aware mode switching. Libraries such as RBDL or Pinocchio can be used to compute \( \mathbf{M}(q) \) and \( \mathbf{h}(q,v) \).

#include <iostream>
#include <vector>

enum class Mode { Free, Contact };

struct State {
    std::vector<double> q;
    std::vector<double> v;
};

struct ContactInfo {
    bool active;
    double gap;        // signed distance
    double normalVel;  // end-effector normal velocity
};

class ImpactController {
public:
    ImpactController(double e) : e_(e), mode_(Mode::Free) {}

    void step(double t, double dt, State &x) {
        ContactInfo c = detectContact(x);
        if (mode_ == Mode::Free && c.gap <= 0.0 && c.normalVel < 0.0) {
            // impact event
            applyImpact(x, c);
            mode_ = Mode::Contact;
        } else if (mode_ == Mode::Contact && c.normalVel >= 0.0 && c.active == false) {
            // liftoff
            mode_ = Mode::Free;
        }

        std::vector<double> tau = computeControl(x, mode_);
        integrateDynamics(t, dt, x, tau);
    }

private:
    Mode mode_;
    double e_; // restitution

    ContactInfo detectContact(const State &x) {
        ContactInfo c{};
        // user-defined kinematics to compute gap and normal velocity
        // c.gap = phi(q); c.normalVel = dphi(q, v);
        return c;
    }

    void applyImpact(State &x, const ContactInfo &c) {
        // Use impact law: v_plus = v_minus - M^{-1} J^T * Lambda_n
        // Here we only sketch the structure; actual implementation requires
        // access to M(q) and J_n(q).
        // M * (v_plus - v_minus) = J_n^T * Lambda_n
        // dot_xn_minus = J_n * v_minus
        // Lambda_n = -(1+e) * dot_xn_minus / (J_n * M^{-1} * J_n^T)
    }

    std::vector<double> computeControl(const State &x, Mode mode) {
        // PD or model-based control with mode-dependent damping gains
        std::vector<double> tau(x.q.size(), 0.0);
        // ...
        return tau;
    }

    void integrateDynamics(double t, double dt, State &x,
                           const std::vector<double> &tau) {
        // Use your favourite integrator (Euler, Runge-Kutta, etc.)
        // M(q) * dv/dt + h(q, v) = tau + J_c^T f_c
    }
};
      

6.2 Java Skeleton (Hybrid Mode Handling)

In Java, we may encode the hybrid modes as an enum and define a controller class that handles transitions. Linear algebra can be handled via libraries like EJML, but here we keep the code scalar for the 1-DOF case.

public class OneDofImpactController {

    public enum Mode { FREE, CONTACT }

    private final double m;
    private final double g;
    private final double e;
    private Mode mode;

    public OneDofImpactController(double m, double g, double e) {
        this.m = m;
        this.g = g;
        this.e = e;
        this.mode = Mode.FREE;
    }

    public static class State {
        public double q;
        public double v;
    }

    public void step(double t, double dt, State x, double qDes) {
        // detect events
        if (mode == Mode.FREE && x.q <= 0.0 && x.v < 0.0) {
            // impact
            x.v = -e * x.v;
            x.q = 0.0;
            mode = Mode.CONTACT;
        } else if (mode == Mode.CONTACT && x.v >= 0.0 && x.q >= 0.0) {
            // liftoff
            mode = Mode.FREE;
        }

        double tau;
        if (mode == Mode.FREE) {
            tau = pdControl(x.q, x.v, qDes, 0.0);
        } else {
            // contact mode: hold at q = 0 with high damping
            tau = 0.0;
        }

        // integrate with explicit Euler
        double a = (tau - m * g) / m;
        x.q = x.q + dt * x.v;
        x.v = x.v + dt * a;
    }

    private double pdControl(double q, double v, double qDes, double vDes) {
        double kp = 50.0;
        double kd = 10.0;
        double tauG = m * g;
        double tauPD = -kp * (q - qDes) - kd * (v - vDes);
        return tauG + tauPD;
    }
}
      

7.1 MATLAB Script with Event Handling

MATLAB supports event functions to detect contact times during numerical integration. Below is a 1-DOF example similar to the Python system.

function contact_transition_demo
    m = 1.0; g = 9.81; e = 0.2;
    x0 = [0.5; -1.0];  % initial [q; v]
    tspan = [0 3];

    opts = odeset('Events', @(t, x) impactEvent(t, x));
    [t, x] = ode45(@(t, x) dynamics(t, x, m, g), tspan, x0, opts);

    % If an impact occurs, apply reset and continue (simple loop)
    tAll = t; xAll = x;
    while t(end) < tspan(2)
        xImpact = x(end, :).';
        if xImpact(1) == 0 && xImpact(2) < 0
            xImpact(2) = -e * xImpact(2); % velocity reset
        end
        tspan2 = [t(end) tspan(2)];
        [t, x, te, xe, ie] = ode45(@(t, x) dynamics(t, x, m, g), ...
                                   tspan2, xImpact, opts);
        tAll = [tAll; t(2:end)];
        xAll = [xAll; x(2:end, :)];
        if isempty(te)
            break;
        end
    end

    figure;
    subplot(2, 1, 1); plot(tAll, xAll(:, 1)); ylabel('q');
    subplot(2, 1, 2); plot(tAll, xAll(:, 2)); ylabel('v'); xlabel('t');

end

function dx = dynamics(t, x, m, g)
    q = x(1); v = x(2);
    tau = pd_control(q, v, 0.2, 0.0, m, g);
    a = (tau - m * g) / m;
    dx = [v; a];
end

function [value, isterminal, direction] = impactEvent(t, x)
    q = x(1); v = x(2);
    value = q;          % detect q = 0
    isterminal = 1;     % stop integration
    direction = -1;     % only when crossing downward
end

function tau = pd_control(q, v, q_des, v_des, m, g)
    kp = 50; kd = 10;
    tau_g = m * g;
    tau_pd = -kp * (q - q_des) - kd * (v - v_des);
    tau = tau_g + tau_pd;
end
      

In Simulink, this logic can be implemented using an ODE solver block, a Hit Crossing block to detect \( q = 0 \), and a triggered subsystem that applies the velocity reset.

7.2 Wolfram Mathematica: Symbolic Impact Map

Mathematica allows closed-form derivation of the impact law and direct numerical simulation using NDSolve with WhenEvent.

(* Define symbols *)
Clear[m, g, e, q, v, t];
m = 1.0; g = 9.81; e = 0.3;

(* Dynamics for free motion *)
eqs = {
  q'[t] == v[t],
  v'[t] == -g
};

(* Impact event: q(t) == 0 and v(t) < 0 *)
impactRule = {v[t] :> -e v[t]};

sol = NDSolve[
  {
    eqs,
    q[0] == 0.5,
    v[0] == -1.0,
    WhenEvent[q[t] == 0 && v[t] < 0,
      {v[t] -> -e v[t], q[t] -> 0}
    ]
  },
  {q, v},
  {t, 0, 3}
];

Plot[{q[t] /. sol[[1]]}, {t, 0, 3},
  PlotLegends -> {"q(t)"}
]
      

For multidimensional manipulators, the same pattern extends to vector states with impact rules that implement the matrix formula derived in Section 2.

In this lesson we modeled contact transitions and impacts as hybrid phenomena in robot control. Starting from impulse-momentum conservation, we derived the rigid-body impact law with a coefficient of restitution, and we showed that kinetic energy cannot increase during impact. We embedded these ideas into a hybrid system formalism with guard sets and reset maps, and we discussed impact-aware control strategies such as pre-impact velocity shaping and post-impact damping modulation. Finally, we illustrated implementation patterns across Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica for simple 1-DOF systems, laying the groundwork for contact-safe controllers in more complex manipulators.