In Chapter 4 Lesson 1, we established that dynamics matter whenever actuator limits, contact forces, or acceleration transients affect motion. Slip and skid are precisely contact-force phenomena: the wheel–ground interface saturates, deforms, and dissipates energy, so the commanded wheel motion does not uniquely determine the body motion.

Let the wheel frame define the longitudinal axis \( x \) and lateral axis \( y \). Denote the wheel angular speed by \( \omega \), wheel radius by \( R \), and contact-frame longitudinal/lateral ground speed by \( v_x, v_y \). Under ideal rolling (no slip), the tangential velocity match is \( v_x = R\omega \) and \( v_y = 0 \), which produces the classic rolling constraints from kinematics. Slip/skid arise when these equalities no longer hold.

flowchart TD
  A["Given: wheel type, speed regime, terrain"] --> B["Define slip measures (kappa, alpha)"]
  B --> C["Choose interaction law"]
  C --> C1["Hard ground: friction-limited \n(Coulomb/brush + saturation)"]
  C --> C2["Soft terrain: \nsinkage + shear deformation \n(pressure-sinkage, shear law)"]
  C --> D["Compute contact forces \nFx, Fy, Fz"]
  D --> E["Integrate dynamics: \nm*dv/dt, I*domega/dt"]
  E --> F["Validate: friction bounds, \nenergy dissipation, plausibility"]
        

The modeling strategy above is the core workflow used in AMR simulation, controller robustness analysis, and design trade studies. The remainder of this lesson makes each block explicit and mathematically defensible.

Longitudinal slip (traction/braking). A common symmetric definition of slip ratio is

\[ \kappa \;=\; \frac{R\omega - v_x}{\max\left(|v_x|,\;|R\omega|,\;\varepsilon\right)}, \quad \varepsilon > 0. \]

Here \( \kappa > 0 \) indicates driving slip (wheel tends to spin faster than ground speed), while \( \kappa < 0 \) indicates braking slip. The max-denominator prevents numerical blow-up near zero speed and preserves scale invariance.

Lateral slip (skid). For small angles, lateral behavior is often captured by a slip angle

\[ \alpha \;=\; \arctan\!\left(\frac{v_y}{|v_x|+\varepsilon}\right). \]

In a skid-steer turn, lateral slip can be substantial even if \( \omega \) tracks commands, because wheel contact forces saturate and the platform “scrubs” laterally.

Proposition (pure rolling equivalence). Assume \( R > 0 \) and \( \max(|v_x|,|R\omega|) > 0 \). Then \( \kappa = 0 \) if and only if \( v_x = R\omega \).

Proof. Under the stated assumption, the denominator in the definition of \( \kappa \) is strictly positive. Therefore \( \kappa = 0 \iff R\omega - v_x = 0 \iff v_x = R\omega \). Conversely, if \( v_x = R\omega \) then the numerator is zero, so \( \kappa=0 \). ■

This proposition is the formal statement that the rolling constraint is recovered as a special case of the slip model.

On rigid terrain, the key physical constraint is that tangential contact force magnitude is bounded by friction. Let the normal load be \( F_z \ge 0 \) and friction coefficient be \( \mu \ge 0 \). A minimal physically consistent bound is \( \sqrt{F_x^2 + F_y^2} \le \mu F_z \).

3.1 Coulomb friction (velocity-direction model). For a relative slip velocity \( \mathbf{v}_{rel} = [v_{rel,x}, v_{rel,y}]^\top \), Coulomb friction can be written as

\[ \mathbf{F}_t \;=\; -\mu F_z \frac{\mathbf{v}_{rel}}{\|\mathbf{v}_{rel}\|+\varepsilon}. \]

This is non-smooth as \( \|\mathbf{v}_{rel}\|\rightarrow 0 \) and often too “stiff” for control-oriented simulation. Therefore, many AMR stacks use regularized laws.

3.2 Linear slip + saturation (control-oriented brush approximation). In the small-slip regime, tire/roller contact can be approximated by linear relations

\[ F_x^{lin} = C_{\kappa}\kappa,\qquad F_y^{lin} = -C_{\alpha}\alpha, \]

with stiffnesses \( C_{\kappa}, C_{\alpha} > 0 \). Physical consistency is enforced by saturating the result so that the friction limit holds.

A widely used hard-ground constraint is the friction ellipse (or circle in normalized coordinates):

\[ \mathcal{E} \;=\; \left\{(F_x,F_y)\;:\; \left(\frac{F_x}{\mu F_z}\right)^2 + \left(\frac{F_y}{\mu F_z}\right)^2 \le 1 \right\}. \]

A common implementation is: compute a “desired” tangential force from a linear model, then project it onto \( \mathcal{E} \).

Theorem (projection onto friction ellipse equals radial scaling). Let \( \mathbf{F}^{lin} = [F_x^{lin},F_y^{lin}]^\top \) and define \( \mathbf{u} = \left[\frac{F_x^{lin}}{\mu F_z}, \frac{F_y^{lin}}{\mu F_z}\right]^\top \), with \( \mu F_z > 0 \). Consider the optimization problem

\[ \min_{\mathbf{F}\in\mathbb{R}^2}\;\|\mathbf{F}-\mathbf{F}^{lin}\|^2 \quad \text{s.t.}\quad \left(\frac{F_x}{\mu F_z}\right)^2 + \left(\frac{F_y}{\mu F_z}\right)^2 \le 1. \]

If \( \|\mathbf{u}\| \le 1 \) then \( \mathbf{F}^\star=\mathbf{F}^{lin} \). Otherwise

\[ \mathbf{F}^\star \;=\; \frac{1}{\|\mathbf{u}\|}\,\mathbf{F}^{lin} \;=\; \frac{\mu F_z}{\sqrt{(F_x^{lin})^2+(F_y^{lin})^2}}\,\mathbf{F}^{lin}. \]

Proof (KKT conditions). Define normalized decision \( \mathbf{w} = \left[\frac{F_x}{\mu F_z}, \frac{F_y}{\mu F_z}\right]^\top \) and note that \( \mathbf{F} = \mu F_z \mathbf{w} \). The objective becomes \( \|\mu F_z \mathbf{w} - \mu F_z \mathbf{u}\|^2 = (\mu F_z)^2 \|\mathbf{w}-\mathbf{u}\|^2 \), so the problem is equivalent to projecting \( \mathbf{u} \) onto the unit ball \( \{\mathbf{w}:\|\mathbf{w}\|\le 1\} \).

The Lagrangian is

\[ \mathcal{L}(\mathbf{w},\lambda)=\|\mathbf{w}-\mathbf{u}\|^2 + \lambda(\|\mathbf{w}\|^2-1), \quad \lambda \ge 0. \]

Stationarity gives \( 2(\mathbf{w}-\mathbf{u}) + 2\lambda \mathbf{w}=0 \Rightarrow (1+\lambda)\mathbf{w}=\mathbf{u} \). If \( \|\mathbf{u}\|\le 1 \), the unconstrained minimizer \( \mathbf{w}=\mathbf{u} \) satisfies feasibility, so \( \lambda=0 \) and \( \mathbf{w}^\star=\mathbf{u} \).

If \( \|\mathbf{u}\| > 1 \), complementary slackness implies the constraint is active: \( \|\mathbf{w}^\star\|=1 \). Using stationarity, \( \mathbf{w}^\star = \mathbf{u}/(1+\lambda) \), and enforcing \( \|\mathbf{w}^\star\|=1 \) yields \( 1+\lambda=\|\mathbf{u}\| \), hence \( \mathbf{w}^\star = \mathbf{u}/\|\mathbf{u}\| \). Returning to forces gives the stated radial scaling. ■

This theorem justifies a simple and numerically robust saturation method: compute linear forces, then apply a single scalar scale factor whenever the friction ellipse would be violated.

On deformable terrain, forces depend not only on instantaneous velocity but also on soil compaction and shear deformation. Full terramechanics models integrate stresses over the wheel–soil contact patch. Here we introduce the core constitutive laws used for those integrals.

5.1 Pressure–sinkage (Bekker-type) relation. For sinkage depth \( z \ge 0 \) and effective contact width \( b > 0 \), a classical relation is

\[ p(z) \;=\; \left(\frac{k_c}{b} + k_{\phi}\right) z^n, \quad k_c \ge 0,\; k_{\phi}\ge 0,\; n > 0, \]

where \( p \) is normal pressure. The normal load is obtained by integrating \( p \) over the contact area:

\[ F_z \;=\; \int_{A} p\!\left(z(\xi)\right)\, dA. \]

5.2 Shear stress vs shear displacement (Janosi–Hanamoto). Let \( c \) be soil cohesion, \( \phi \) internal friction angle, and \( \sigma(\xi) \) normal stress at patch coordinate \( \xi \). The shear stress magnitude can be modeled as

\[ \tau(\xi) \;=\; \left(c + \sigma(\xi)\tan\phi\right)\left(1-\exp\!\left(-\frac{j(\xi)}{K}\right)\right), \quad K > 0, \]

where \( j(\xi) \) is the local shear displacement and \( K \) is a shear deformation modulus. The limiting cases are informative:

\[ \begin{aligned} j(\xi) \approx 0 \;⇒\; \tau(\xi) \approx \left(c+\sigma(\xi)\tan\phi\right)\frac{j(\xi)}{K} \quad \text{(approximately linear)}\\[4pt] j(\xi)\rightarrow\infty \;⇒\; \tau(\xi)\rightarrow c+\sigma(\xi)\tan\phi \quad \text{(saturation)} \end{aligned} \]

5.3 A tractable surrogate for AMR simulation. Full integration requires a geometric contact model, \( z(\xi) \), and a shear displacement law \( j(\xi) \) as a function of slip. In many control-level AMR simulations, a compact surrogate is used that preserves the essential behavior: near-linear traction at small slip and saturation at high slip:

\[ F_x \;=\; F_z\,\mu_{eff}(\kappa),\qquad \mu_{eff}(\kappa) \;=\; \mu_{peak}\left(1-e^{-k|\kappa|}\right)\operatorname{sgn}(\kappa), \quad \mu_{peak}\in(0,1),\; k>0. \]

This surrogate is not a replacement for terramechanics integrals, but it is widely used to make slip-aware controllers and planners testable without high computational cost.

The practical value of slip models depends on parameter identification. We outline two common identification tasks that do not require concepts beyond standard least squares.

6.1 Hard-ground linear stiffness (small-slip data). Suppose we can measure \( (F_x,\kappa) \) pairs in a region where saturation is not active: \( |F_x| \ll \mu F_z \). Then \( F_x \approx C_{\kappa}\kappa \) and

\[ \hat{C}_{\kappa} \;=\; \arg\min_{C_{\kappa}} \sum_{i=1}^{N}\left(F_{x,i}-C_{\kappa}\kappa_i\right)^2 \;=\; \frac{\sum_{i=1}^{N} \kappa_i F_{x,i}}{\sum_{i=1}^{N}\kappa_i^2}, \quad \text{if } \sum \kappa_i^2 > 0. \]

6.2 Pressure–sinkage regression (plate/wheel tests). If we have \( (z_i,p_i) \) pairs and assume \( p_i = K z_i^n \) with \( K=\left(\frac{k_c}{b}+k_{\phi}\right) \), then taking logs gives

\[ \log p_i \;=\; \log K + n\log z_i, \]

which reduces to linear regression in \( (\log z_i,\log p_i) \). This is the standard starting point for terramechanics parameter estimation.

flowchart TD
  A["Collect data: encoders, IMU, load estimate"] --> B["Compute slip measures: kappa, alpha"]
  B --> C["Filter to small-slip region (no saturation)"]
  C --> D["Fit Ck, Ca by least squares"]
  B --> E["If soft terrain: measure (z,p) from tests"]
  E --> F["Fit K and n via log-linear regression"]
  D --> G["Validate: predict Fx,Fy; check friction bounds"]
  F --> G
        

File: Chapter4_Lesson2.py

"""
Chapter4_Lesson2.py
Autonomous Mobile Robots — Chapter 4 Lesson 2
Slip, Skid, and Terrain Interaction Models

This script provides:
1) Hard-ground wheel slip model with combined-slip saturation (friction ellipse).
2) Simple soft-soil (terramechanics-inspired) traction surrogate.
3) A small time-domain simulation of a driven wheel + 1D vehicle mass.

Dependencies: numpy, matplotlib (standard scientific Python stack).
"""

from __future__ import annotations
import numpy as np
import matplotlib.pyplot as plt


def slip_ratio(vx: float, omega: float, R: float, eps: float = 1e-6) -> float:
    """
    Longitudinal slip ratio kappa in a symmetric definition:

        kappa = (R*omega - vx) / max(|vx|, |R*omega|, eps)

    Positive kappa corresponds to driving (wheel tends to spin faster than ground speed).
    Negative kappa corresponds to braking (wheel tends to rotate slower than ground speed).
    """
    denom = max(abs(vx), abs(R * omega), eps)
    return (R * omega - vx) / denom


def slip_angle(vx: float, vy: float, eps: float = 1e-6) -> float:
    """
    Lateral slip angle alpha = atan2(vy, |vx| + eps).
    """
    return np.arctan2(vy, abs(vx) + eps)


def friction_ellipse_saturate(Fx_lin: float, Fy_lin: float, mu: float, Fz: float) -> tuple[float, float]:
    """
    Enforce (Fx/(mu Fz))^2 + (Fy/(mu Fz))^2 <= 1 by radial scaling.

    This is the Euclidean projection onto a scaled L2 ball in normalized coordinates,
    which equals radial scaling for this particular constraint set.
    """
    cap = max(mu * Fz, 0.0)
    if cap <= 0.0:
        return 0.0, 0.0

    nx = Fx_lin / cap
    ny = Fy_lin / cap
    r2 = nx * nx + ny * ny
    if r2 <= 1.0:
        return Fx_lin, Fy_lin

    s = 1.0 / np.sqrt(r2)
    return s * Fx_lin, s * Fy_lin


def hard_ground_forces(vx: float, vy: float, omega: float, R: float,
                       Fz: float, mu: float,
                       Ck: float = 15000.0, Ca: float = 12000.0) -> tuple[float, float, float, float]:
    """
    Simple hard-ground tire model:
      Fx = Ck * kappa,  Fy = -Ca * alpha  (small-slip linearization)
      then saturate with friction ellipse.

    Returns: Fx, Fy, kappa, alpha
    """
    kappa = slip_ratio(vx, omega, R)
    alpha = slip_angle(vx, vy)
    Fx_lin = Ck * kappa
    Fy_lin = -Ca * alpha
    Fx, Fy = friction_ellipse_saturate(Fx_lin, Fy_lin, mu, Fz)
    return Fx, Fy, kappa, alpha


def soft_soil_traction_surrogate(kappa: float, Fz: float, mu_peak: float = 0.65, k_shape: float = 8.0) -> float:
    """
    A compact terramechanics-inspired surrogate for longitudinal traction:
        Fx = Fz * mu_eff(kappa),
        mu_eff(kappa) = mu_peak * (1 - exp(-k_shape * |kappa|)) * sign(kappa)

    This captures:
      - approximately linear traction near kappa = 0,
      - saturation to a peak traction coefficient as |kappa| increases.

    It is NOT a full Bekker + Janosi-Hanamoto integral model, but is useful for control-level simulation.
    """
    s = np.sign(kappa) if kappa != 0.0 else 0.0
    mu_eff = mu_peak * (1.0 - np.exp(-k_shape * abs(kappa))) * s
    return Fz * mu_eff


def simulate_1d_wheel_vehicle(
    T_cmd: float = 12.0,
    terrain: str = "hard",
    t_end: float = 4.0,
    dt: float = 1e-3,
) -> dict[str, np.ndarray]:
    """
    1D model: vehicle longitudinal speed vx, wheel spin omega.

    Dynamics:
      m * dvx/dt = Fx - Frr
      Iw * domega/dt = T_cmd - R*Fx - bw*omega

    For 'hard' terrain: Fx from linear tire + friction ellipse.
    For 'soft' terrain: Fx from traction surrogate.

    Returns time histories in a dictionary.
    """
    # Parameters (typical small AMR)
    m = 25.0           # kg
    R = 0.10           # m
    Iw = 0.05          # kg m^2 (wheel+gear equivalent)
    bw = 0.02          # N m s/rad viscous loss
    g = 9.81

    # Normal load per wheel (assume one driven wheel supporting part of mass)
    Fz = 0.25 * m * g

    # Rolling resistance
    Crr = 0.02
    Frr = Crr * Fz

    # Hard-ground friction
    mu = 0.8

    N = int(t_end / dt) + 1
    t = np.linspace(0.0, t_end, N)
    vx = np.zeros(N)
    omega = np.zeros(N)
    kappa = np.zeros(N)
    Fx = np.zeros(N)

    # Initial conditions
    vx[0] = 0.0
    omega[0] = 0.0

    for i in range(N - 1):
        if terrain == "hard":
            Fx_i, _, k_i, _ = hard_ground_forces(vx[i], 0.0, omega[i], R, Fz, mu)
        elif terrain == "soft":
            k_i = slip_ratio(vx[i], omega[i], R)
            Fx_i = soft_soil_traction_surrogate(k_i, Fz, mu_peak=0.55, k_shape=10.0)
            # Soil also increases motion resistance: raise effective rolling resistance a bit
            Fx_i -= 0.015 * Fz * np.sign(vx[i])
        else:
            raise ValueError("terrain must be 'hard' or 'soft'.")

        # Vehicle
        ax = (Fx_i - Frr) / m
        vx[i + 1] = vx[i] + dt * ax

        # Wheel
        domega = (T_cmd - R * Fx_i - bw * omega[i]) / Iw
        omega[i + 1] = omega[i] + dt * domega

        Fx[i] = Fx_i
        kappa[i] = k_i

    Fx[-1] = Fx[-2]
    kappa[-1] = kappa[-2]

    return dict(t=t, vx=vx, omega=omega, kappa=kappa, Fx=Fx, R=np.array([R]), Fz=np.array([Fz]))


def main() -> None:
    hard = simulate_1d_wheel_vehicle(T_cmd=12.0, terrain="hard")
    soft = simulate_1d_wheel_vehicle(T_cmd=12.0, terrain="soft")

    # Plot
    plt.figure()
    plt.plot(hard["t"], hard["vx"], label="vx hard")
    plt.plot(soft["t"], soft["vx"], label="vx soft")
    plt.xlabel("time [s]")
    plt.ylabel("vx [m/s]")
    plt.grid(True)
    plt.legend()

    plt.figure()
    plt.plot(hard["t"], hard["kappa"], label="kappa hard")
    plt.plot(soft["t"], soft["kappa"], label="kappa soft")
    plt.xlabel("time [s]")
    plt.ylabel("slip ratio kappa [-]")
    plt.grid(True)
    plt.legend()

    plt.figure()
    plt.plot(hard["t"], hard["Fx"], label="Fx hard")
    plt.plot(soft["t"], soft["Fx"], label="Fx soft")
    plt.xlabel("time [s]")
    plt.ylabel("Fx [N]")
    plt.grid(True)
    plt.legend()

    plt.show()


if __name__ == "__main__":
    main()
      

File: Chapter4_Lesson2.cpp

/*
Chapter4_Lesson2.cpp
Autonomous Mobile Robots — Chapter 4 Lesson 2
Slip, Skid, and Terrain Interaction Models

This C++ example mirrors the Python script:
- Hard-ground combined-slip saturation with a friction ellipse
- Simple 1D wheel + vehicle simulation

Suggested libraries in robotics stacks:
  - Eigen (linear algebra) for larger models
  - ROS 2 (rclcpp) for real-time integration (not used here)

Build (example):
  g++ -O2 -std=c++17 Chapter4_Lesson2.cpp -o lesson2

*/

#include <algorithm>
#include <cmath>
#include <iostream>
#include <string>
#include <vector>

static double slip_ratio(double vx, double omega, double R, double eps = 1e-6) {
    const double denom = std::max({std::abs(vx), std::abs(R * omega), eps});
    return (R * omega - vx) / denom;
}

static double friction_ellipse_scale(double Fx, double Fy, double mu, double Fz) {
    const double cap = std::max(mu * Fz, 0.0);
    if (cap <= 0.0) return 0.0;
    const double nx = Fx / cap;
    const double ny = Fy / cap;
    const double r2 = nx * nx + ny * ny;
    if (r2 <= 1.0) return 1.0;
    return 1.0 / std::sqrt(r2);
}

static double hard_ground_Fx(double vx, double omega, double R, double Fz, double mu,
                             double Ck = 15000.0) {
    const double kappa = slip_ratio(vx, omega, R);
    const double Fx_lin = Ck * kappa;
    // 1D: ignore lateral force, so ellipse reduces to |Fx| <= mu Fz
    const double cap = std::max(mu * Fz, 0.0);
    return std::clamp(Fx_lin, -cap, cap);
}

static double soft_soil_Fx(double kappa, double Fz, double mu_peak = 0.55, double k_shape = 10.0) {
    const double s = (kappa > 0.0) ? 1.0 : (kappa < 0.0 ? -1.0 : 0.0);
    const double mu_eff = mu_peak * (1.0 - std::exp(-k_shape * std::abs(kappa))) * s;
    return Fz * mu_eff;
}

struct SimResult {
    std::vector<double> t, vx, omega, kappa, Fx;
};

static SimResult simulate_1d(double T_cmd, const std::string& terrain, double t_end = 4.0, double dt = 1e-3) {
    // Parameters
    const double m = 25.0;
    const double R = 0.10;
    const double Iw = 0.05;
    const double bw = 0.02;
    const double g = 9.81;
    const double Fz = 0.25 * m * g;
    const double Crr = 0.02;
    const double Frr = Crr * Fz;
    const double mu = 0.8;

    const int N = static_cast<int>(t_end / dt) + 1;

    SimResult out;
    out.t.resize(N);
    out.vx.assign(N, 0.0);
    out.omega.assign(N, 0.0);
    out.kappa.assign(N, 0.0);
    out.Fx.assign(N, 0.0);

    for (int i = 0; i < N; ++i) out.t[i] = i * dt;

    for (int i = 0; i < N - 1; ++i) {
        double Fx_i = 0.0;
        double k_i = slip_ratio(out.vx[i], out.omega[i], R);

        if (terrain == "hard") {
            Fx_i = hard_ground_Fx(out.vx[i], out.omega[i], R, Fz, mu);
        } else if (terrain == "soft") {
            Fx_i = soft_soil_Fx(k_i, Fz);
            Fx_i -= 0.015 * Fz * ((out.vx[i] > 0.0) ? 1.0 : (out.vx[i] < 0.0 ? -1.0 : 0.0));
        } else {
            throw std::runtime_error("terrain must be 'hard' or 'soft'");
        }

        const double ax = (Fx_i - Frr) / m;
        out.vx[i + 1] = out.vx[i] + dt * ax;

        const double domega = (T_cmd - R * Fx_i - bw * out.omega[i]) / Iw;
        out.omega[i + 1] = out.omega[i] + dt * domega;

        out.Fx[i] = Fx_i;
        out.kappa[i] = k_i;
    }

    out.Fx[N - 1] = out.Fx[N - 2];
    out.kappa[N - 1] = out.kappa[N - 2];
    return out;
}

int main() {
    try {
        auto hard = simulate_1d(12.0, "hard");
        auto soft = simulate_1d(12.0, "soft");

        // Print a small summary (final values)
        const auto idx = static_cast<int>(hard.t.size()) - 1;
        std::cout << "Final (hard): vx=" << hard.vx[idx] << " m/s, kappa=" << hard.kappa[idx] << ", Fx=" << hard.Fx[idx] << " N\n";
        std::cout << "Final (soft): vx=" << soft.vx[idx] << " m/s, kappa=" << soft.kappa[idx] << ", Fx=" << soft.Fx[idx] << " N\n";

        std::cout << "Tip: export vectors to CSV if you want to plot externally.\n";
    } catch (const std::exception& e) {
        std::cerr << "Error: " << e.what() << "\n";
        return 1;
    }
    return 0;
}
      

File: Chapter4_Lesson2.java

/*
Chapter4_Lesson2.java
Autonomous Mobile Robots — Chapter 4 Lesson 2
Slip, Skid, and Terrain Interaction Models

This Java example provides a minimal 1D simulation (vehicle vx, wheel omega).
In Java robotics ecosystems, common libraries include:
  - EJML for linear algebra (not required here)
  - ROS 2 Java client libraries for integration (not used here)

Compile and run:
  javac Chapter4_Lesson2.java
  java Chapter4_Lesson2
*/

import java.util.Arrays;

public class Chapter4_Lesson2 {

    static double slipRatio(double vx, double omega, double R, double eps) {
        double denom = Math.max(Math.max(Math.abs(vx), Math.abs(R * omega)), eps);
        return (R * omega - vx) / denom;
    }

    static double hardGroundFx(double vx, double omega, double R, double Fz, double mu, double Ck) {
        double kappa = slipRatio(vx, omega, R, 1e-6);
        double FxLin = Ck * kappa;
        double cap = Math.max(mu * Fz, 0.0);
        return Math.max(-cap, Math.min(cap, FxLin));
    }

    static double softSoilFx(double kappa, double Fz, double muPeak, double kShape) {
        double s = (kappa > 0.0) ? 1.0 : (kappa < 0.0 ? -1.0 : 0.0);
        double muEff = muPeak * (1.0 - Math.exp(-kShape * Math.abs(kappa))) * s;
        return Fz * muEff;
    }

    static class SimResult {
        double[] t, vx, omega, kappa, Fx;
        SimResult(int N) {
            t = new double[N];
            vx = new double[N];
            omega = new double[N];
            kappa = new double[N];
            Fx = new double[N];
        }
    }

    static SimResult simulate1D(double Tcmd, String terrain, double tEnd, double dt) {
        // Parameters
        double m = 25.0;
        double R = 0.10;
        double Iw = 0.05;
        double bw = 0.02;
        double g = 9.81;

        double Fz = 0.25 * m * g;
        double Crr = 0.02;
        double Frr = Crr * Fz;
        double mu = 0.8;

        int N = (int)Math.round(tEnd / dt) + 1;
        SimResult out = new SimResult(N);

        for (int i = 0; i < N; i++) out.t[i] = i * dt;

        for (int i = 0; i < N - 1; i++) {
            double k = slipRatio(out.vx[i], out.omega[i], R, 1e-6);
            double Fx;

            if ("hard".equals(terrain)) {
                Fx = hardGroundFx(out.vx[i], out.omega[i], R, Fz, mu, 15000.0);
            } else if ("soft".equals(terrain)) {
                Fx = softSoilFx(k, Fz, 0.55, 10.0);
                Fx -= 0.015 * Fz * Math.signum(out.vx[i]);
            } else {
                throw new IllegalArgumentException("terrain must be 'hard' or 'soft'");
            }

            double ax = (Fx - Frr) / m;
            out.vx[i + 1] = out.vx[i] + dt * ax;

            double domega = (Tcmd - R * Fx - bw * out.omega[i]) / Iw;
            out.omega[i + 1] = out.omega[i] + dt * domega;

            out.Fx[i] = Fx;
            out.kappa[i] = k;
        }

        out.Fx[N - 1] = out.Fx[N - 2];
        out.kappa[N - 1] = out.kappa[N - 2];
        return out;
    }

    public static void main(String[] args) {
        SimResult hard = simulate1D(12.0, "hard", 4.0, 1e-3);
        SimResult soft = simulate1D(12.0, "soft", 4.0, 1e-3);

        int idx = hard.t.length - 1;
        System.out.println("Final (hard): vx=" + hard.vx[idx] + " m/s, kappa=" + hard.kappa[idx] + ", Fx=" + hard.Fx[idx] + " N");
        System.out.println("Final (soft): vx=" + soft.vx[idx] + " m/s, kappa=" + soft.kappa[idx] + ", Fx=" + soft.Fx[idx] + " N");

        System.out.println("If you want plots, write arrays to CSV and plot in your preferred tool.");
    }
}
      

File: Chapter4_Lesson2.m

% Chapter4_Lesson2.m
% Autonomous Mobile Robots — Chapter 4 Lesson 2
% Slip, Skid, and Terrain Interaction Models
%
% This MATLAB script includes:
%  (1) hard-ground slip model (linear + saturation),
%  (2) a terramechanics-inspired traction surrogate,
%  (3) an OPTIONAL programmatic Simulink model builder for the same 1D dynamics.
%
% Requirements:
%   - MATLAB
%   - Simulink (only for the Simulink build section)

clear; clc;

%% Parameters
m  = 25.0;      % kg
R  = 0.10;      % m
Iw = 0.05;      % kg*m^2
bw = 0.02;      % N*m*s/rad
g  = 9.81;
Fz = 0.25*m*g;

mu  = 0.8;      % hard-ground friction
Ck  = 15000.0;  % longitudinal stiffness
Crr = 0.02;     % rolling resistance coefficient
Frr = Crr*Fz;

%% Simulation settings
Tcmd = 12.0;    % N*m
dt = 1e-3;
tEnd = 4.0;
t = 0:dt:tEnd;

[vx_h, om_h, kap_h, Fx_h] = simulate_1d(t, dt, Tcmd, "hard", m, R, Iw, bw, Fz, mu, Ck, Frr);
[vx_s, om_s, kap_s, Fx_s] = simulate_1d(t, dt, Tcmd, "soft", m, R, Iw, bw, Fz, mu, Ck, Frr);

%% Plots
figure; plot(t, vx_h, t, vx_s); grid on;
xlabel('time [s]'); ylabel('vx [m/s]'); legend('hard','soft');

figure; plot(t, kap_h, t, kap_s); grid on;
xlabel('time [s]'); ylabel('slip ratio kappa [-]'); legend('hard','soft');

figure; plot(t, Fx_h, t, Fx_s); grid on;
xlabel('time [s]'); ylabel('Fx [N]'); legend('hard','soft');

%% OPTIONAL: Build a Simulink model programmatically
% This section creates a Simulink model "Chapter4_Lesson2_Sim.slx" that simulates
% the same 1D dynamics using Integrator blocks and a MATLAB Function block.
%
% Comment out if you do not have Simulink.
try
    build_simulink_model();
catch ME
    fprintf('Simulink model build skipped or failed: %s\n', ME.message);
end

%% -------- Local functions --------

function kappa = slip_ratio(vx, omega, R)
    eps = 1e-6;
    denom = max([abs(vx), abs(R*omega), eps]);
    kappa = (R*omega - vx)/denom;
end

function Fx = hard_ground_Fx(vx, omega, R, Fz, mu, Ck)
    kappa = slip_ratio(vx, omega, R);
    Fx_lin = Ck*kappa;
    cap = max(mu*Fz, 0);
    Fx = min(max(Fx_lin, -cap), cap);
end

function Fx = soft_soil_Fx(kappa, Fz)
    muPeak = 0.55; kShape = 10.0;
    s = sign(kappa);
    muEff = muPeak*(1 - exp(-kShape*abs(kappa)))*s;
    Fx = Fz*muEff;
end

function [vx, omega, kappa_hist, Fx_hist] = simulate_1d(t, dt, Tcmd, terrain, m, R, Iw, bw, Fz, mu, Ck, Frr)
    N = length(t);
    vx = zeros(1,N);
    omega = zeros(1,N);
    kappa_hist = zeros(1,N);
    Fx_hist = zeros(1,N);

    for i=1:N-1
        k = slip_ratio(vx(i), omega(i), R);
        if terrain == "hard"
            Fx = hard_ground_Fx(vx(i), omega(i), R, Fz, mu, Ck);
        else
            Fx = soft_soil_Fx(k, Fz);
            Fx = Fx - 0.015*Fz*sign(vx(i)); % extra soil motion resistance surrogate
        end

        ax = (Fx - Frr)/m;
        vx(i+1) = vx(i) + dt*ax;

        domega = (Tcmd - R*Fx - bw*omega(i))/Iw;
        omega(i+1) = omega(i) + dt*domega;

        Fx_hist(i) = Fx;
        kappa_hist(i) = k;
    end
    Fx_hist(end) = Fx_hist(end-1);
    kappa_hist(end) = kappa_hist(end-1);
end

function build_simulink_model()
    model = 'Chapter4_Lesson2_Sim';
    if bdIsLoaded(model)
        close_system(model, 0);
    end
    new_system(model); open_system(model);

    % Add blocks
    add_block('simulink/Sources/Constant', [model '/Tcmd'], 'Value', '12');
    add_block('simulink/Continuous/Integrator', [model '/Int_vx'], 'InitialCondition', '0');
    add_block('simulink/Continuous/Integrator', [model '/Int_omega'], 'InitialCondition', '0');

    add_block('simulink/User-Defined Functions/MATLAB Function', [model '/ForcesAndDerivs']);
    add_block('simulink/Sinks/Scope', [model '/Scope']);

    % Parameters as workspace variables (use base workspace)
    assignin('base','m',25.0);
    assignin('base','R',0.10);
    assignin('base','Iw',0.05);
    assignin('base','bw',0.02);
    assignin('base','g',9.81);
    assignin('base','Fz',0.25*25.0*9.81);
    assignin('base','mu',0.8);
    assignin('base','Ck',15000.0);
    assignin('base','Crr',0.02);
    assignin('base','Frr',0.02*(0.25*25.0*9.81));
    assignin('base','terrainFlag',1); % 1=hard, 2=soft

    % Configure MATLAB Function block code
    code = [
        "function [dvx, domega, kappa, Fx] = f(vx, omega, Tcmd)", newline, ...
        "%#codegen", newline, ...
        "eps = 1e-6;", newline, ...
        "den = max([abs(vx), abs(R*omega), eps]);", newline, ...
        "kappa = (R*omega - vx)/den;", newline, ...
        "if terrainFlag == 1", newline, ...
        "  Fx_lin = Ck*kappa;", newline, ...
        "  cap = max(mu*Fz,0);", newline, ...
        "  Fx = min(max(Fx_lin,-cap),cap);", newline, ...
        "else", newline, ...
        "  muPeak = 0.55; kShape = 10.0;", newline, ...
        "  Fx = Fz*(muPeak*(1-exp(-kShape*abs(kappa)))*sign(kappa));", newline, ...
        "  Fx = Fx - 0.015*Fz*sign(vx);", newline, ...
        "end", newline, ...
        "dvx = (Fx - Frr)/m;", newline, ...
        "domega = (Tcmd - R*Fx - bw*omega)/Iw;", newline ...
    ];
    set_param([model '/ForcesAndDerivs'], 'Script', code);

    % Wire blocks
    % Tcmd -> ForcesAndDerivs
    add_line(model, 'Tcmd/1', 'ForcesAndDerivs/3');
    % vx -> ForcesAndDerivs
    add_line(model, 'Int_vx/1', 'ForcesAndDerivs/1');
    % omega -> ForcesAndDerivs
    add_line(model, 'Int_omega/1', 'ForcesAndDerivs/2');

    % dvx -> Int_vx
    add_line(model, 'ForcesAndDerivs/1', 'Int_vx/1');
    % domega -> Int_omega
    add_line(model, 'ForcesAndDerivs/2', 'Int_omega/1');

    % outputs to Scope: vx, omega, kappa, Fx
    add_block('simulink/Signal Routing/Mux', [model '/Mux'], 'Inputs', '4');
    add_line(model, 'Int_vx/1', 'Mux/1');
    add_line(model, 'Int_omega/1', 'Mux/2');
    add_line(model, 'ForcesAndDerivs/3', 'Mux/3');
    add_line(model, 'ForcesAndDerivs/4', 'Mux/4');
    add_line(model, 'Mux/1', 'Scope/1');

    % Layout (rough)
    set_param([model '/Tcmd'], 'Position', [50 80 100 110]);
    set_param([model '/Int_vx'], 'Position', [280 40 310 70]);
    set_param([model '/Int_omega'], 'Position', [280 120 310 150]);
    set_param([model '/ForcesAndDerivs'], 'Position', [160 60 250 140]);
    set_param([model '/Mux'], 'Position', [360 55 390 155]);
    set_param([model '/Scope'], 'Position', [430 80 460 130]);

    % Simulation settings
    set_param(model, 'StopTime', '4');
    save_system(model);
    fprintf('Created Simulink model: %s.slx\n', model);

    % Run once
    sim(model);
end
      

File: Chapter4_Lesson2.nb

(* ::Package:: *)

(* Chapter4_Lesson2.nb
   Autonomous Mobile Robots — Chapter 4 Lesson 2
   Slip, Skid, and Terrain Interaction Models

   This notebook-style file (plain text) can be opened in Wolfram Mathematica.
*)

Notebook[{
  Cell["Chapter 4 — Lesson 2: Slip, Skid, and Terrain Interaction Models", "Title"],

  Cell["1) Slip ratio and basic properties", "Section"],
  Cell[BoxData @ ToBoxes @ HoldForm[
    kappa[vx_, omega_, R_, eps_: 10^-6] := (R*omega - vx)/Max[Abs[vx], Abs[R*omega], eps]
  ], "Input"],

  Cell["Show: kappa = 0 corresponds to pure rolling (vx = R omega when both nonzero).", "Text"],
  Cell[BoxData @ ToBoxes @ HoldForm[
    Simplify[kappa[vx, omega, R] == 0, {vx != 0, omega != 0, R > 0}]
  ], "Input"],

  Cell["2) Hard-ground force model with friction ellipse", "Section"],
  Cell[BoxData @ ToBoxes @ HoldForm[
    saturateEllipse[{Fx_, Fy_}, mu_, Fz_] := Module[{cap = mu*Fz, r2},
      If[cap <= 0, {0, 0},
        r2 = (Fx/cap)^2 + (Fy/cap)^2;
        If[r2 <= 1, {Fx, Fy}, {Fx, Fy}/Sqrt[r2]]
      ]
    ]
  ], "Input"],

  Cell["3) Small simulation (1D wheel + vehicle) using NDSolve", "Section"],
  Cell[BoxData @ ToBoxes @ HoldForm[
    params = <|"m" -> 25, "R" -> 0.10, "Iw" -> 0.05, "bw" -> 0.02, "g" -> 9.81, "mu" -> 0.8, 
              "Ck" -> 15000, "Crr" -> 0.02, "Tcmd" -> 12|>;
    Fz = 0.25*params["m"]*params["g"];
    Frr = params["Crr"]*Fz;

    FxHard[vx_, omega_] := Module[{kap, FxLin, cap},
      kap = kappa[vx, omega, params["R"]];
      FxLin = params["Ck"]*kap;
      cap = params["mu"]*Fz;
      Clip[FxLin, {-cap, cap}]
    ];

    sol = NDSolve[{
        vx'[t] == (FxHard[vx[t], omega[t]] - Frr)/params["m"],
        omega'[t] == (params["Tcmd"] - params["R"]*FxHard[vx[t], omega[t]] - params["bw"]*omega[t])/params["Iw"],
        vx[0] == 0, omega[0] == 0
      },
      {vx, omega}, {t, 0, 4}
    ][[1]];

    Plot[Evaluate[{vx[t], omega[t]} /. sol], {t, 0, 4}, PlotLegends -> {"vx(t)", "omega(t)"}]
  ], "Input"]
}]
      

Problem 1 (Slip ratio boundedness): Using the definition \( \kappa = \frac{R\omega - v_x}{\max(|v_x|,|R\omega|,\varepsilon)} \), show that \( |\kappa|\le 2 \) for any \( v_x,\omega \). Under what conditions does \( |\kappa|\le 1 \) hold?

Solution: Let \( d=\max(|v_x|,|R\omega|,\varepsilon) \). Then \( |\kappa| = \frac{|R\omega - v_x|}{d} \le \frac{|R\omega|+|v_x|}{d} \). Since \( d \ge |R\omega| \) and \( d \ge |v_x| \), we have

\[ |\kappa| \;\le\; \frac{|R\omega|}{d} + \frac{|v_x|}{d} \;\le\; 1 + 1 \;=\; 2. \]

If additionally one uses the traction-only convention with \( |v_x| \ge \varepsilon \) and defines \( \kappa = \frac{R\omega - v_x}{|v_x|} \) for driving (or the braking variant), then \( |\kappa|\le 1 \) holds for common operating ranges where \( R\omega \) and \( v_x \) remain comparable. The symmetric definition guarantees numerical robustness but allows values up to 2 in extreme zero-speed mismatches.

Problem 2 (Friction ellipse feasibility): Prove that if \( (F_x,F_y)\in \mathcal{E} \) then the tangential power dissipation satisfies \( |P_t| \le \mu F_z \|\mathbf{v}_{rel}\| \) for Coulomb friction, where \( P_t = \mathbf{F}_t^\top \mathbf{v}_{rel} \).

Solution: Under Coulomb friction, \( \mathbf{F}_t = -\mu F_z \frac{\mathbf{v}_{rel}}{\|\mathbf{v}_{rel}\|+\varepsilon} \). Then

\[ P_t = \mathbf{F}_t^\top \mathbf{v}_{rel} = -\mu F_z \frac{\mathbf{v}_{rel}^\top \mathbf{v}_{rel}}{\|\mathbf{v}_{rel}\|+\varepsilon} = -\mu F_z \frac{\|\mathbf{v}_{rel}\|^2}{\|\mathbf{v}_{rel}\|+\varepsilon}. \]

Hence \( P_t \le 0 \) (pure dissipation), and \( |P_t| = \mu F_z \frac{\|\mathbf{v}_{rel}\|^2}{\|\mathbf{v}_{rel}\|+\varepsilon} \le \mu F_z \|\mathbf{v}_{rel}\| \). The ellipse constraint provides the same bound via Cauchy–Schwarz: \( |P_t|\le \|\mathbf{F}_t\|\|\mathbf{v}_{rel}\| \le \mu F_z \|\mathbf{v}_{rel}\| \).

Problem 3 (Estimating \(C_{\kappa}\) from data): Given measured pairs \( (\kappa_i,F_{x,i}) \) in the non-saturated region, derive the least squares estimator for \( C_{\kappa} \) and show it is unbiased if \( F_{x,i} = C_{\kappa}\kappa_i + \eta_i \) with \( E[\eta_i]=0 \) and \( \kappa_i \) deterministic.

Solution: Minimizing \( \sum (F_{x,i}-C_{\kappa}\kappa_i)^2 \) and differentiating w.r.t. \( C_{\kappa} \) yields

\[ \hat{C}_{\kappa} = \frac{\sum \kappa_i F_{x,i}}{\sum \kappa_i^2}. \]

Substitute the measurement model: \( \hat{C}_{\kappa} = \frac{\sum \kappa_i(C_{\kappa}\kappa_i+\eta_i)}{\sum \kappa_i^2} = C_{\kappa} + \frac{\sum \kappa_i\eta_i}{\sum \kappa_i^2} \). Taking expectation gives \( E[\hat{C}_{\kappa}] = C_{\kappa} + \frac{\sum \kappa_i E[\eta_i]}{\sum \kappa_i^2} = C_{\kappa} \), so the estimator is unbiased.

Problem 4 (Soft-soil shear law limits): Using \( \tau = (c+\sigma\tan\phi)\left(1-e^{-j/K}\right) \), show the small-shear approximation and the saturation limit.

Solution: For small \( j \), use the Taylor expansion \( e^{-j/K} \approx 1 - j/K \), hence \( 1-e^{-j/K} \approx j/K \). Therefore

\[ \tau \approx (c+\sigma\tan\phi)\frac{j}{K}. \]

For \( j\rightarrow\infty \), \( e^{-j/K}\rightarrow 0 \), so \( \tau \rightarrow c+\sigma\tan\phi \), i.e., shear stress saturates at a Mohr–Coulomb-like limit.

Problem 5 (Steady-state slip in a 1D wheel–vehicle model): Consider the simplified dynamics \( m\dot v_x = F_x - F_{rr} \) and \( I_w \dot\omega = T - R F_x - b_w\omega \). Suppose hard-ground traction saturates at \( F_x = \mu F_z \) in steady driving. Derive the steady-state wheel speed \( \omega^\star \) and vehicle acceleration \( \dot v_x^\star \).

Solution: If saturation holds, then in steady-state wheel speed (constant \( \omega^\star \)) we set \( \dot\omega = 0 \):

\[ 0 = T - R(\mu F_z) - b_w \omega^\star \quad \Rightarrow \quad \omega^\star = \frac{T - R\mu F_z}{b_w}. \]

The vehicle acceleration is then

\[ \dot v_x^\star = \frac{\mu F_z - F_{rr}}{m}. \]

This highlights a key AMR design point: once traction saturates, increasing torque mainly increases wheel spin (hence slip) rather than forward acceleration.

We defined longitudinal slip ratio and lateral slip angle as the minimal quantities that relax ideal rolling constraints. For hard ground, we built a linear slip-force approximation and proved that enforcing the friction ellipse can be done by a principled Euclidean projection (radial scaling). For soft terrain, we introduced pressure–sinkage and shear-deformation constitutive laws and provided a tractable traction surrogate suitable for AMR simulation. These models will be used in subsequent lessons to construct simple full-robot dynamic models and to analyze how parameters affect real navigation behavior.