Odometry for wheeled robots is a deterministic integration of encoder increments under rolling constraints. In practice, the mapping from encoder ticks to ground displacement is biased by effective wheel radii, unequal radii (left vs right), and uncertainty in wheel separation. These biases create systematic drift that cannot be removed by filtering alone without correcting the underlying kinematic parameters.

We model the effective parameters (including tire deformation and consistent slip regimes) as \( r_L \), \( r_R \), and \( b \). The calibration goal is to infer these parameters from recorded motions where we can measure the robot’s relative motion in the plane.

flowchart TD
  A["Encoder increments (dphi_L, dphi_R) + external relative pose (dx, dy, dtheta)"] --> B["Choose parameterization p = [rL, rR, b]"]
  B --> C["Predict body-frame increment from kinematics"]
  C --> D["Residual e(p) = predicted - measured"]
  D --> E["Nonlinear least squares: minimize sum ||e||^2"]
  E --> F["Solve with Gauss-Newton / Levenberg-Marquardt"]
  F --> G["Validate on new trajectories + update nominal parameters"]
        

From Lesson 2 (Differential Drive Kinematics), for a single motion segment (or sample interval) we assume we know the wheel rotation increments \( \Delta\phi_L \) and \( \Delta\phi_R \). The corresponding arc-length increments are \( \Delta s_L \) and \( \Delta s_R \), modeled as:

\[ \Delta s_L = r_L \Delta\phi_L,\qquad \Delta s_R = r_R \Delta\phi_R. \]

Define the average forward increment and yaw increment: \( \Delta s \) and \( \Delta\theta \).

\[ \Delta s = \frac{\Delta s_R + \Delta s_L}{2},\qquad \Delta\theta = \frac{\Delta s_R - \Delta s_L}{b}. \]

For calibration we use a body-frame relative motion model (increment expressed in the robot frame at the start of the segment). Under constant curvature during the segment, the exact body-frame increment is:

\[ \Delta x = \Delta s \,\mathrm{sinc}(\Delta\theta),\qquad \Delta y = \Delta s \,\mathrm{cosc}(\Delta\theta),\qquad \Delta\theta = \frac{r_R\Delta\phi_R - r_L\Delta\phi_L}{b}, \]

where the two scalar functions are defined as:

\[ \mathrm{sinc}(u) = \frac{\sin(u)}{u},\qquad \mathrm{cosc}(u) = \frac{1-\cos(u)}{u}, \]

with continuous extensions at \( u=0 \) by limits:

\[ \lim_{u \rightarrow 0}\mathrm{sinc}(u)=1,\qquad \lim_{u \rightarrow 0}\mathrm{cosc}(u)=0. \]

Proof (limits): Using Taylor expansions \( \sin(u)=u-\frac{u^3}{6}+\mathcal{O}(u^5) \) and \( \cos(u)=1-\frac{u^2}{2}+\mathcal{O}(u^4) \), we obtain:

\[ \mathrm{sinc}(u)=\frac{u-\frac{u^3}{6}+\mathcal{O}(u^5)}{u} =1-\frac{u^2}{6}+\mathcal{O}(u^4), \qquad \mathrm{cosc}(u)=\frac{\frac{u^2}{2}+\mathcal{O}(u^4)}{u} =\frac{u}{2}+\mathcal{O}(u^3). \]

Hence \( \Delta x \) and \( \Delta y \) are numerically stable if we evaluate these functions with small-angle safeguards (implemented in the code sections).

flowchart LR
  WL["Left wheel: rL, dphi_L"] --- C["Chassis center"]
  C --- WR["Right wheel: rR, dphi_R"]
  WL --- B["baseline b"] --- WR
  C --> T["Body increment (dx, dy, dtheta) from kinematics"]
        

Let \( p = [r_L, r_R, b]^\top \). Each segment provides three measured quantities \( (\Delta x^{gt},\Delta y^{gt},\Delta\theta^{gt}) \) (from an external source such as motion capture, AprilTag tracking, or a tape-measured maneuver), and two encoder inputs \( (\Delta\phi_L,\Delta\phi_R) \).

Key point: If the dataset does not excite rotation (i.e., \( \Delta\theta \approx 0 \) for all segments), then \( b \) is not identifiable because it only appears in \( \Delta\theta = \frac{\Delta s_R-\Delta s_L}{b} \).

A useful “minimal” pair of maneuvers illustrates the structure:

(i) Straight-line maneuver. If ground truth indicates \( \Delta\theta^{gt}=0 \), then the model implies:

\[ 0 = \Delta\theta = \frac{r_R\Delta\phi_R - r_L\Delta\phi_L}{b} \quad \Longrightarrow \quad \frac{r_R}{r_L} = \frac{\Delta\phi_L}{\Delta\phi_R}. \]

This identifies only the ratio \( r_R/r_L \), not the absolute scale.

(ii) In-place rotation maneuver. If the robot rotates about its center with negligible translation, then \( \Delta s_R \approx -\Delta s_L \), and ground truth provides \( \Delta\theta^{gt} \neq 0 \). From

\[ \Delta\theta^{gt} \approx \frac{\Delta s_R - \Delta s_L}{b} = \frac{r_R\Delta\phi_R - r_L\Delta\phi_L}{b}, \]

once \( r_R/r_L \) is known (from straight motion), the remaining scale and \( b \) can be determined from sufficiently many rotation segments (or one well-measured rotation with known wheel increments).

Rank condition (local identifiability): For a dataset of \( N \) segments, stack residuals into \( e(p)\in\mathbb{R}^{3N} \). If the Jacobian \( J(p)=\frac{\partial e}{\partial p}\in\mathbb{R}^{3N\times 3} \) has full column rank (rank 3) near the solution, then \( p \) is locally identifiable and Gauss–Newton has a well-defined normal equation.

For segment \( k \), define the predicted increment: \( \hat{\Delta x}_k(p),\hat{\Delta y}_k(p),\hat{\Delta\theta}_k(p) \) from the model in Section 2, and measured ground-truth increment: \( \Delta x_k^{gt},\Delta y_k^{gt},\Delta\theta_k^{gt} \).

Define the residual vector per segment:

\[ e_k(p) = \begin{bmatrix} \hat{\Delta x}_k(p) - \Delta x_k^{gt}\\ \hat{\Delta y}_k(p) - \Delta y_k^{gt}\\ \mathrm{wrap}\!\left(\hat{\Delta\theta}_k(p) - \Delta\theta_k^{gt}\right) \end{bmatrix}. \]

The calibration objective (with optional weights \( W_k \succ 0 \)) is:

\[ p^\star = \arg\min_{p} \sum_{k=1}^{N} e_k(p)^\top W_k e_k(p). \]

In the common unweighted case \( W_k = I_3 \), this is standard nonlinear least squares. If translation is measured more accurately than orientation (or vice versa), weights improve conditioning.

Gauss–Newton step: Linearize around \( p \): \( e(p+\delta)\approx e(p)+J(p)\delta \). The Gauss–Newton update solves:

\[ \delta^\star = \arg\min_{\delta}\|e(p)+J(p)\delta\|_2^2 \quad \Longrightarrow \quad (J^\top J)\delta^\star = -J^\top e. \]

Levenberg–Marquardt damping: To handle mild nonlinearity and improve robustness, solve:

\[ (J^\top J + \lambda I)\delta^\star = -J^\top e,\qquad \lambda > 0. \]

Proof (descent direction when \( J^\top e \neq 0 \)): Let \( \delta = -(J^\top J + \lambda I)^{-1}J^\top e \). Since \( J^\top J + \lambda I \succ 0 \), we have \( \delta^\top J^\top e = - e^\top J (J^\top J + \lambda I)^{-1} J^\top e \le 0 \), with strict inequality if \( J^\top e \neq 0 \). Thus for sufficiently small step size, the objective decreases.

Analytic Jacobians significantly improve convergence and numerical stability. Define:

\[ s_L = r_L\Delta\phi_L,\quad s_R = r_R\Delta\phi_R,\quad \Delta s = \frac{s_R+s_L}{2},\quad \Delta\theta = \frac{s_R-s_L}{b}. \]

Using \( \Delta x = \Delta s\,\mathrm{sinc}(\Delta\theta) \) and \( \Delta y = \Delta s\,\mathrm{cosc}(\Delta\theta) \), we apply the chain rule. First compute partial derivatives of \( \Delta s \) and \( \Delta\theta \).

\[ \frac{\partial \Delta s}{\partial r_L} = \frac{\Delta\phi_L}{2},\qquad \frac{\partial \Delta s}{\partial r_R} = \frac{\Delta\phi_R}{2},\qquad \frac{\partial \Delta s}{\partial b} = 0, \]

\[ \frac{\partial \Delta\theta}{\partial r_L} = -\frac{\Delta\phi_L}{b},\qquad \frac{\partial \Delta\theta}{\partial r_R} = \frac{\Delta\phi_R}{b},\qquad \frac{\partial \Delta\theta}{\partial b} = -\frac{s_R-s_L}{b^2}=-\frac{\Delta\theta}{b}. \]

Next compute derivatives of the special functions. Let \( f(u)=\mathrm{sinc}(u)=\frac{\sin u}{u} \), \( g(u)=\mathrm{cosc}(u)=\frac{1-\cos u}{u} \). For \( u \neq 0 \):

\[ f'(u) = \frac{u\cos(u)-\sin(u)}{u^2},\qquad g'(u) = \frac{u\sin(u)-(1-\cos(u))}{u^2}. \]

Proof (quotient rule): For \( f(u)=\frac{\sin u}{u} \), \( f'(u)=\frac{u\cos u - \sin u}{u^2} \). For \( g(u)=\frac{1-\cos u}{u} \), \( g'(u)=\frac{u\sin u -(1-\cos u)}{u^2} \). At \( u=0 \), use Taylor expansions to define stable series approximations.

Finally:

\[ \frac{\partial \Delta x}{\partial p_i} = \frac{\partial \Delta s}{\partial p_i} f(\Delta\theta) + \Delta s\, f'(\Delta\theta)\frac{\partial \Delta\theta}{\partial p_i}, \qquad \frac{\partial \Delta y}{\partial p_i} = \frac{\partial \Delta s}{\partial p_i} g(\Delta\theta) + \Delta s\, g'(\Delta\theta)\frac{\partial \Delta\theta}{\partial p_i}. \]

These expressions are implemented directly in all provided language implementations.

A robust calibration dataset should include a mixture of straight segments, arcs with varying curvature, and in-place rotations. This guarantees Jacobian rank and improves conditioning.

Minimum recommendations:

• At least 30–50 segments with nonzero \( |\Delta\theta| \) (turning excitation).
• At least 30 segments with near-straight motion (to stabilize scale of \( \Delta s \)).
• Ground-truth increments measured in the robot body frame at each segment start (or converted to it).

Common pitfalls:

• Unmodeled slip: if slip varies strongly across the dataset, the effective parameters are not constant; the optimizer will “average” and may not generalize.
• Poor angle wrapping: always compare orientation increments using a wrapping function to keep residuals in \( (-\pi,\pi] \).
• Unit mismatch: ensure encoder increments are in radians and translations in meters.

This implementation uses numpy and (optionally) scipy.optimize.least_squares. It also includes a toolbox-free Gauss–Newton fallback. The dataset format is a CSV with columns: dphi_L, dphi_R, dx_gt, dy_gt, dtheta_gt.

Chapter2_Lesson5.py

#!/usr/bin/env python3
# Chapter2_Lesson5.py
"""
Chapter 2 — Wheeled Locomotion Kinematics (Mobile-Specific)
Lesson 5 — Kinematic Calibration of Wheel Parameters

This script performs nonlinear least-squares calibration of differential-drive wheel parameters:
  p = [r_L, r_R, b]
from paired data:
  (dphi_L[k], dphi_R[k]) and ground-truth relative increments (dx_gt[k], dy_gt[k], dtheta_gt[k])
where (dx, dy, dtheta) are expressed in the robot body frame at the beginning of each segment.

Dependencies (recommended):
  - numpy
  - scipy  (scipy.optimize.least_squares)
"""

from __future__ import annotations

import csv
from dataclasses import dataclass
from typing import Tuple

import numpy as np

try:
    from scipy.optimize import least_squares
except Exception:
    least_squares = None


@dataclass
class CalibData:
    dphi_L: np.ndarray
    dphi_R: np.ndarray
    dx_gt: np.ndarray
    dy_gt: np.ndarray
    dtheta_gt: np.ndarray


def wrap_to_pi(a: np.ndarray) -> np.ndarray:
    return (a + np.pi) % (2.0 * np.pi) - np.pi


def _sinc1(x: np.ndarray) -> np.ndarray:
    y = np.empty_like(x)
    small = np.abs(x) < 1e-6
    xs = x[small]
    y[small] = 1.0 - (xs * xs) / 6.0 + (xs ** 4) / 120.0
    xb = x[~small]
    y[~small] = np.sin(xb) / xb
    return y


def _cosc1(x: np.ndarray) -> np.ndarray:
    y = np.empty_like(x)
    small = np.abs(x) < 1e-6
    xs = x[small]
    y[small] = xs / 2.0 - (xs ** 3) / 24.0 + (xs ** 5) / 720.0
    xb = x[~small]
    y[~small] = (1.0 - np.cos(xb)) / xb
    return y


def _dsinc1(x: np.ndarray) -> np.ndarray:
    y = np.empty_like(x)
    small = np.abs(x) < 1e-5
    xs = x[small]
    y[small] = -(xs) / 3.0 + (xs ** 3) / 30.0 - (xs ** 5) / 840.0
    xb = x[~small]
    y[~small] = (xb * np.cos(xb) - np.sin(xb)) / (xb * xb)
    return y


def _dcosc1(x: np.ndarray) -> np.ndarray:
    y = np.empty_like(x)
    small = np.abs(x) < 1e-5
    xs = x[small]
    y[small] = 0.5 - (xs * xs) / 8.0 + (xs ** 4) / 144.0
    xb = x[~small]
    y[~small] = (xb * np.sin(xb) - (1.0 - np.cos(xb))) / (xb * xb)
    return y


def predict_body_increment(dphi_L: np.ndarray, dphi_R: np.ndarray, p: np.ndarray) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    rL, rR, b = float(p[0]), float(p[1]), float(p[2])
    sL = rL * dphi_L
    sR = rR * dphi_R
    ds = 0.5 * (sR + sL)
    dtheta = (sR - sL) / b
    dx = ds * _sinc1(dtheta)
    dy = ds * _cosc1(dtheta)
    return dx, dy, dtheta


def residuals_and_jacobian(p: np.ndarray, data: CalibData) -> Tuple[np.ndarray, np.ndarray]:
    rL, rR, b = float(p[0]), float(p[1]), float(p[2])
    dphiL = data.dphi_L
    dphiR = data.dphi_R

    sL = rL * dphiL
    sR = rR * dphiR
    A = sR + sL
    B = sR - sL

    ds = 0.5 * A
    dtheta = B / b

    f = _sinc1(dtheta)
    g = _cosc1(dtheta)
    fp = _dsinc1(dtheta)
    gp = _dcosc1(dtheta)

    dx = ds * f
    dy = ds * g

    ex = dx - data.dx_gt
    ey = dy - data.dy_gt
    eth = wrap_to_pi(dtheta - data.dtheta_gt)

    dds_drL = 0.5 * dphiL
    dds_drR = 0.5 * dphiR

    ddth_drL = -(dphiL) / b
    ddth_drR = (dphiR) / b
    ddth_db = -(B) / (b * b)

    ddx_drL = dds_drL * f + ds * fp * ddth_drL
    ddx_drR = dds_drR * f + ds * fp * ddth_drR
    ddx_db = 0.0 * f + ds * fp * ddth_db

    ddy_drL = dds_drL * g + ds * gp * ddth_drL
    ddy_drR = dds_drR * g + ds * gp * ddth_drR
    ddy_db = 0.0 * g + ds * gp * ddth_db

    N = dphiL.size
    r = np.empty(3 * N, dtype=float)
    J = np.empty((3 * N, 3), dtype=float)

    r[0::3] = ex
    r[1::3] = ey
    r[2::3] = eth

    J[0::3, 0] = ddx_drL
    J[0::3, 1] = ddx_drR
    J[0::3, 2] = ddx_db

    J[1::3, 0] = ddy_drL
    J[1::3, 1] = ddy_drR
    J[1::3, 2] = ddy_db

    J[2::3, 0] = ddth_drL
    J[2::3, 1] = ddth_drR
    J[2::3, 2] = ddth_db

    return r, J


def load_csv(path: str) -> CalibData:
    rows = []
    with open(path, "r", newline="") as f:
        reader = csv.reader(f)
        for row in reader:
            if not row:
                continue
            if row[0].strip().lower().startswith("dphi"):
                continue
            rows.append([float(x) for x in row[:5]])
    arr = np.asarray(rows, dtype=float)
    return CalibData(
        dphi_L=arr[:, 0],
        dphi_R=arr[:, 1],
        dx_gt=arr[:, 2],
        dy_gt=arr[:, 3],
        dtheta_gt=arr[:, 4],
    )


def calibrate(data: CalibData, p0: np.ndarray) -> np.ndarray:
    if least_squares is not None:
        def fun(p):
            r, _ = residuals_and_jacobian(p, data)
            return r

        def jac(p):
            _, J = residuals_and_jacobian(p, data)
            return J

        res = least_squares(fun, x0=p0, jac=jac, method="trf", loss="huber", f_scale=0.01)
        return res.x.astype(float)

    # Fallback: damped Gauss-Newton
    p = p0.astype(float).copy()
    lam = 1e-3
    for _ in range(50):
        r, J = residuals_and_jacobian(p, data)
        A = J.T @ J + lam * np.eye(3)
        g = J.T @ r
        dp = -np.linalg.solve(A, g)
        if np.linalg.norm(dp) < 1e-10:
            break
        p_new = p + dp
        r_new, _ = residuals_and_jacobian(p_new, data)
        if (r_new @ r_new) < (r @ r):
            p = p_new
            lam *= 0.7
        else:
            lam *= 2.0
    return p


def main() -> None:
    import argparse
    parser = argparse.ArgumentParser()
    parser.add_argument("--csv", type=str, default="")
    parser.add_argument("--rL0", type=float, default=0.05)
    parser.add_argument("--rR0", type=float, default=0.05)
    parser.add_argument("--b0", type=float, default=0.30)
    args = parser.parse_args()

    if args.csv:
        data = load_csv(args.csv)
    else:
        # minimal synthetic example
        rng = np.random.default_rng(7)
        N = 200
        dphiL = rng.normal(0.0, 0.4, N)
        dphiR = rng.normal(0.0, 0.4, N)
        p_true = np.array([0.050, 0.051, 0.300], dtype=float)
        dx, dy, dth = predict_body_increment(dphiL, dphiR, p_true)
        data = CalibData(dphiL, dphiR, dx + 0.002*rng.normal(size=N), dy + 0.002*rng.normal(size=N), dth + 0.001*rng.normal(size=N))

    p0 = np.array([args.rL0, args.rR0, args.b0], dtype=float)
    phat = calibrate(data, p0=p0)

    print("Estimated [rL, rR, b] =", phat)

if __name__ == "__main__":
    main()
      

Robotics software note: In ROS/ROS2, the calibrated parameters typically update the wheel radius and wheel separation used by the odometry node; you then regenerate nav_msgs/Odometry (or equivalent) from encoder streams using the corrected kinematics.

The C++ implementation uses Eigen (commonly used throughout robotics stacks) to solve the 3-parameter normal equations. For production systems, the same residual and Jacobian can be used with Ceres Solver; the core math remains unchanged.

Chapter2_Lesson5.cpp

// Chapter2_Lesson5.cpp
#include <Eigen/Dense>
#include <cmath>
#include <fstream>
#include <iostream>
#include <sstream>
#include <string>
#include <vector>

/* (Full source is provided in the downloadable ZIP.) */
      

The downloadable source includes: CSV parsing, stable sinc/cosc evaluations, analytic Jacobians, and a Levenberg-style damping strategy.

Java is common in industrial tooling and some academic environments. We use EJML for matrix operations, mirroring the Eigen-based approach.

Chapter2_Lesson5.java

// Chapter2_Lesson5.java
// (Full source is provided in the downloadable ZIP.)
// It includes stable sinc/cosc functions, analytic Jacobians, CSV parsing,
// and a damped Gauss-Newton solver using EJML SimpleMatrix.
      

MATLAB is widely used for controls and robotics prototyping. The provided script implements a toolbox-free damped Gauss–Newton solver and includes an optional function to programmatically build a small Simulink model computing \( (\Delta x,\Delta y,\Delta\theta) \) from wheel increments and parameters.

Chapter2_Lesson5.m

% Chapter2_Lesson5.m
% (Full source is provided in the downloadable ZIP.)
% Includes: simulateDataset, residualAndJacobian, calibrateGN, and buildSimulinkModel.
      

Mathematica is valuable for symbolic manipulation and rapid verification. The provided notebook expression defines the same residual model and runs a robust quasi-Newton minimization.

Chapter2_Lesson5.nb

(* Chapter2_Lesson5.nb *)
(* (Full notebook expression is provided in the downloadable ZIP.) *)
      

Problem 1 (Straight-Line Constraint and Radius Ratio): A differential-drive robot performs a straight-line motion segment for which ground truth indicates \( \Delta\theta^{gt}=0 \). Show that the wheel radii satisfy a ratio constraint depending only on \( \Delta\phi_L \) and \( \Delta\phi_R \), and explain why the absolute scale is not identifiable from straight segments alone.

Solution: From the kinematic model:

\[ 0=\Delta\theta=\frac{r_R\Delta\phi_R-r_L\Delta\phi_L}{b} \quad \Longrightarrow \quad r_R\Delta\phi_R=r_L\Delta\phi_L \quad \Longrightarrow \quad \frac{r_R}{r_L}=\frac{\Delta\phi_L}{\Delta\phi_R}. \]

This constrains only the ratio \( r_R/r_L \). Multiplying both radii by a constant \( c \) scales \( \Delta s \) and therefore \( \Delta x \) by \( c \), but straight-line data alone provides no independent measurement to fix \( c \) unless translation ground truth is used. If only the constraint \( \Delta\theta^{gt}=0 \) is used (without \( \Delta x^{gt} \)), the scale remains free.

Problem 2 (Local Identifiability via Jacobian Rank): Consider the stacked residual vector \( e(p)\in\mathbb{R}^{3N} \) and Jacobian \( J(p)\in\mathbb{R}^{3N\times 3} \). Provide a concrete condition on the dataset that makes it plausible that \( \mathrm{rank}(J)=3 \). Explain why “only straight segments” typically fails this condition.

Solution: A plausible sufficient condition is that the dataset contains: (i) segments with diverse \( (\Delta\phi_L,\Delta\phi_R) \) producing nonzero \( \Delta\theta \) (turning excitation), and (ii) segments with nonzero average motion \( \Delta s \). If all segments have \( \Delta\theta \approx 0 \), then the residual component for \( \Delta\theta \) becomes nearly insensitive to \( b \) (since \( \Delta\theta=(s_R-s_L)/b \approx 0 \) and \( s_R-s_L \approx 0 \)), making the third column of \( J \) nearly linearly dependent or close to zero across rows. This reduces rank and makes \( b \) poorly identifiable.

Problem 3 (Derive \( f'(u) \) and \( g'(u) \)): Let \( f(u)=\frac{\sin u}{u} \) and \( g(u)=\frac{1-\cos u}{u} \) for \( u\neq 0 \). Derive \( f'(u) \) and \( g'(u) \).

Solution: Use the quotient rule:

\[ f'(u)=\frac{u\cos(u)-\sin(u)}{u^2},\qquad g'(u)=\frac{u\sin(u)-(1-\cos(u))}{u^2}. \]

For numerical stability at \( |u| < \epsilon \), use Taylor-series expansions to evaluate \( f \), \( g \), \( f' \), and \( g' \).

Problem 4 (Curvature Bias from Wheel-Radius Mismatch): Suppose the true wheel radii are \( r_L=r(1-\varepsilon) \) and \( r_R=r(1+\varepsilon) \), while the odometry assumes both are \( r \). For a commanded “straight” segment where encoders report \( \Delta\phi_L=\Delta\phi_R=\Delta\phi \), derive the true yaw increment \( \Delta\theta \) and show it is proportional to \( \varepsilon \).

Solution: The true wheel displacements are:

\[ \Delta s_L = r(1-\varepsilon)\Delta\phi,\qquad \Delta s_R = r(1+\varepsilon)\Delta\phi. \]

Therefore the true yaw increment is:

\[ \Delta\theta = \frac{\Delta s_R - \Delta s_L}{b} = \frac{r(1+\varepsilon)\Delta\phi - r(1-\varepsilon)\Delta\phi}{b} = \frac{2r\varepsilon\Delta\phi}{b}. \]

Hence even when encoder increments are equal, a radius mismatch induces systematic curvature proportional to \( \varepsilon \). Kinematic calibration estimates this mismatch via \( r_L \) and \( r_R \).

We constructed a kinematics-only calibration method for differential-drive wheel parameters by: (i) expressing relative motion increments in the body frame using stable \( \mathrm{sinc} \) and \( \mathrm{cosc} \) forms; (ii) proving key limits and deriving analytic Jacobians; and (iii) formulating calibration as nonlinear least squares solved with Gauss–Newton / Levenberg–Marquardt. The provided implementations (Python, C++, Java, MATLAB/Simulink, Mathematica) directly encode these derivations.