Integrate the squared jerk in closed form rather than sampling it

README.md

@@ -52,6 +52,13 @@ is a comparison quantity, not an exact figure. The minimum-jerk cost has no such
 jerk is continuous: `720 d^2 / T^5 = 720 * 2.25 / 2.34375^5 = 22.9065`, which the test suite checks
 against the closed form directly.
 
+The caveat applies to the table rather than to the library. `certify_jerk_cost` integrates the
+coefficient table segment by segment rather than running a quadrature over a trace, which is exact
+for a piecewise polynomial however discontinuous the integrand is, and for this S-curve it returns
+239.99999999999994, the analytic 240.0 to fifteen digits, at any sample count, because it reads no
+samples at all. The table above keeps the quadrature so that every column of it comes from the same
+instrument.
+
 How much of each budget the three profiles actually consume:
 
 | Profile | Velocity used | Acceleration used | Jerk used |
@@ -96,6 +103,7 @@ The package ships a `py.typed` marker, so it delivers its annotations to anythin
 from trajectory_optimizer import (
     LimitSet,
     PointToPointMove,
+    certify_jerk_cost,
     certify_limits,
     minimum_jerk_time_optimal,
     sample_trajectory,
@@ -109,12 +117,14 @@ report = certify_limits(trajectory, limits)
 
 print(f"duration {trajectory.duration:.5f} s")
 print(f"peak velocity {trajectory.peak_derivatives().row('velocity')[0]:.5f} rad/s")
+print(f"jerk cost {certify_jerk_cost(trajectory):.5f} rad^2/s^5")
 print(f"within limits everywhere: {report.ok}")
 ```
 
 ```text
 duration 2.34375 s
 peak velocity 1.20000 rad/s
+jerk cost 22.90649 rad^2/s^5
 within limits everywhere: True
 ```
 
@@ -288,6 +298,12 @@ instant, and the time scaling search uses it for its feasibility predicate. That
 peaks that matter often sit at irrational fractions of the duration: the acceleration peak of a
 rest-to-rest minimum-jerk move is at `t / T = (1 - 1 / sqrt(3)) / 2`, which no uniform grid contains.
 
+The objective is integrated rather than sampled by the same argument. The square of a segment's jerk
+polynomial is a polynomial too, so the integral over that segment is a quadratic form in the jerk
+coefficients, and `certify_jerk_cost` sums that form over every segment and every joint. Scaling
+divides the result by the fifth power of its factor and synchronisation adds its components, so both
+wrappers forward the exact value rather than resampling what they wrap.
+
 Time-optimal time scaling holds the geometric path fixed and searches for the smallest stretch
 factor whose scaled derivatives all fit inside their bounds. Because a stretch by `s` divides
 velocity by `s`, acceleration by `s` squared, and jerk by `s` cubed, feasibility is monotone in `s`,
@@ -302,8 +318,8 @@ distinction, the alternatives that were rejected, and the limitations that remai
 | Module | Responsibility |
 | --- | --- |
 | `src/trajectory_optimizer/model.py` | Pure dataclasses for joint limits, via points, boundary conditions, and moves, with validation and no input or output |
-| `src/trajectory_optimizer/algorithm/protocol.py` | The `Trajectory` protocol every generator satisfies, the `CertifiedTrajectory` protocol for those that solve for their own peaks, and the single-instant state record |
-| `src/trajectory_optimizer/algorithm/piecewise.py` | Piecewise polynomial container giving exact derivatives and exact peaks, plus the mapping of a scalar path profile onto a joint-space line |
+| `src/trajectory_optimizer/algorithm/protocol.py` | The `Trajectory` protocol every generator satisfies, the `CertifiedTrajectory` and `IntegrableTrajectory` protocols for those that solve for their own peaks and integrate their own jerk cost, and the single-instant state record |
+| `src/trajectory_optimizer/algorithm/piecewise.py` | Piecewise polynomial container giving exact derivatives, exact peaks, and the exact jerk cost integral, plus the mapping of a scalar path profile onto a joint-space line |
 | `src/trajectory_optimizer/algorithm/minimum_jerk.py` | Closed-form quintic minimum-jerk solution and its analytic time-optimal duration |
 | `src/trajectory_optimizer/algorithm/trapezoidal.py` | Trapezoidal velocity profile including the triangular degenerate case |
 | `src/trajectory_optimizer/algorithm/scurve.py` | Seven segment jerk limited profile including both degenerate cases |
@@ -312,7 +328,7 @@ distinction, the alternatives that were rejected, and the limitations that remai
 | `src/trajectory_optimizer/algorithm/synchronisation.py` | Stretching independent joint trajectories onto one common duration |
 | `src/trajectory_optimizer/pipeline.py` | Sampling any trajectory into a structured trace of times and derivatives |
 | `src/trajectory_optimizer/analysis/limits.py` | The sampled limit check, the certified one, and peak usage per joint |
-| `src/trajectory_optimizer/analysis/metrics.py` | Duration, peak derivatives, jerk cost, and joint-space path length |
+| `src/trajectory_optimizer/analysis/metrics.py` | Duration, peak derivatives, joint-space path length, and the jerk cost both as a quadrature of a trace and as an exact integral of the coefficient table |
 | `src/trajectory_optimizer/analysis/report.py` | Markdown rendering of comparison tables, trade-offs, limit reports, and scaling results |
 | `src/trajectory_optimizer/analysis/figures.py` | Stacked derivative figures and profile overlays, drawn without a pyplot backend |
 | `examples/` | Wiring scripts only, with no logic of their own |
@@ -331,8 +347,8 @@ uv run ruff format --check .
 uv run mypy
 ```
 
-The suite is 203 tests over three tiers and completes in about 10 seconds. Coverage of
-`src/trajectory_optimizer` is 97 percent of 1204 statements; CI runs the same command with
+The suite is 218 tests over three tiers and completes in about 10 seconds. Coverage of
+`src/trajectory_optimizer` is 97 percent of 1230 statements; CI runs the same command with
 `--cov-fail-under=95` on Ubuntu and on Windows, along with the linter and mypy in strict mode over
 the package, the examples, and the figure script.
 
@@ -346,6 +362,9 @@ point, and synchronisation gives every joint identical start and end times. The
 checked against `scipy.interpolate.CubicSpline` and the bisection scale factor against the analytic
 one. The closed-form peaks are checked against their analytic expressions, against a two million
 point grid, and against a constructed move whose violation a three point check calls compliant. The
+integrated jerk cost is checked against a hand-computed integral, against `720 d^2 / T^5` for the
+minimum-jerk quintic and `400 * 4 * 0.15` for the S-curve, and against the quadrature it replaces at
+4001, 40001, and 400001 samples, where the error falls by a factor of ten at each step. The
 second tier pins recorded profiles at eleven sampled fractions of each duration to an absolute
 tolerance of 1e-9. The third tier runs every script in `examples/` as a subprocess with
 `--samples 51`.

docs/design-notes.md

@@ -149,6 +149,30 @@ synchronisation concatenates the blocks of its components, which is exact becaus
 owns a disjoint set of joints. Both raise `TypeError` when what they wrap cannot answer, so a
 certificate over an estimate is never produced.
 
+### The jerk cost integrated rather than sampled
+
+The same argument settles the objective itself. On one segment the jerk of one joint is a polynomial
+in the local time, so its square is a polynomial of twice the degree and the integral over that
+segment is a quadratic form in the jerk coefficients `c` of that segment and joint,
+
+```text
+c^T M c,    M_ik = L^(i + k + 1) / (i + k + 1),
+```
+
+with `L` the length of the segment. `PiecewisePolynomialTrajectory.jerk_cost` sums that form over
+every segment and every joint and reads no samples at all, so the value carries no quadrature error
+however discontinuous the integrand is. `certify_jerk_cost` is the analysis-layer entry point.
+
+The two wrappers forward rather than repeat, as they do for the peaks. `TimeScaledTrajectory`
+divides by the fifth power of its scale, because stretching time divides the jerk by the cube of the
+scale and multiplies the interval of integration by it, and `SynchronisedTrajectory` adds its
+components, because the cost is a sum over joints and each component owns a disjoint set of them.
+Both raise `TypeError` when what they wrap cannot answer.
+
+An acceleration that steps at a breakpoint carries an impulse of jerk that the segment polynomials
+do not contain, so the cost reported for a trapezoidal profile is zero. That is the same one-sided
+reading of a piecewise definition as its reported peak jerk of zero, and it is not a bound either.
+
 ### Time-optimal time scaling by bisection
 
 Given a trajectory `q` on `[0, T]`, the scaled trajectory `q_s(t) = q(t / s)` traverses the same
@@ -260,6 +284,36 @@ that is not a piecewise polynomial cannot be certified by this method; its repor
 Third, the certificate covers the kinematic limits only, which is the boundary set out under "No
 dynamics, no actuator model" below.
 
+### The jerk cost was a quadrature over a sampled trace
+
+This entry recorded that the jerk cost was a Simpson quadrature of the sampled squared jerk, which
+converges to machine precision only where the jerk is continuous. For the minimum-jerk profile the
+analytic cost of a 1.5 rad move over 2.0 s is 50.625 and the quadrature returned 50.625000000000014
+at 20001 samples. For the S-curve, whose jerk is piecewise constant and discontinuous, the error
+fell only as the reciprocal of the sample count: against an analytic 240.0 the reported values were
+240.360000, 240.036000, and 240.003600 at 4001, 40001, and 400001 samples, so the figure was a
+comparison quantity rather than an exact one. The entry named its own remedy: integrate the
+coefficient table segment by segment.
+
+That is now implemented and the entry is closed. `PiecewisePolynomialTrajectory.jerk_cost`
+integrates each segment exactly, `certify_jerk_cost` is the entry point in the analysis layer, and
+both wrappers forward. The method is described under "The jerk cost integrated rather than sampled"
+above.
+
+What it cost. Twenty six executable statements across five modules and one new protocol,
+`IntegrableTrajectory`. No new dependency and no new record, because the result is a float. The
+S-curve above now reports 239.99999999999994 at any sample count, which is the analytic 240.0 to
+fifteen digits, and the time-optimal minimum-jerk move reports 22.906492245333652 against a closed
+form of 22.906492245333332. Runtime is one quadratic form of order at most three per segment and
+joint, so it is cheaper than sampling the trace it replaces.
+
+What remains. Two things. First, `compute_metrics` still reads a trace and its jerk cost is still a
+quadrature, because a trajectory that is not a piecewise polynomial cannot be integrated this way;
+its docstring now says so. Second, the path length in the same record is still the polyline through
+the samples, which under-estimates a curved joint-space path by an amount falling as the square of
+the sample spacing. Point to point moves are straight lines in joint space, where the polyline is
+exact, so that error is confined to the spline paths.
+
 ## Known limitations
 
 ### Time scaling along a fixed path is not global time optimality
@@ -308,17 +362,6 @@ segment case analysis does extend to nonzero end velocities, at the cost of seve
 including the case where the profile must first decelerate, and that extension is not implemented.
 Moves with nonzero end derivatives should use the minimum-jerk generator or a clamped spline.
 
-### Metrics are quadratures over sampled traces
-
-The jerk cost is a Simpson quadrature of the sampled squared jerk. For the minimum-jerk profile,
-where jerk is continuous, this converges to machine precision: the analytic cost of a 1.5 rad move
-over 2.0 s is 50.625 and the quadrature returns 50.625000000000014 at 20001 samples. For the
-S-curve, where jerk is piecewise constant and discontinuous, the quadrature error falls only as the
-reciprocal of the sample count: the analytic cost is 240.0 and the reported values are 240.360000,
-240.036000, and 240.003600 at 4001, 40001, and 400001 samples. Reading the S-curve jerk cost as an
-exact figure would be a mistake; it is a comparison quantity. Integrating the coefficient table
-segment by segment would make it exact and is the obvious improvement.
-
 ### No dynamics, no actuator model
 
 Everything here is kinematic. Torque limits, payload-dependent inertia, friction, and joint

src/trajectory_optimizer/__init__.py

@@ -6,6 +6,7 @@
 
 from trajectory_optimizer.algorithm import (
     CertifiedTrajectory,
+    IntegrableTrajectory,
     PeakDerivatives,
     PiecewisePolynomialTrajectory,
     SynchronisedTrajectory,
@@ -21,7 +22,12 @@
     time_optimal_scaling,
     trapezoidal_move,
 )
-from trajectory_optimizer.analysis import certify_limits, check_limits, compute_metrics
+from trajectory_optimizer.analysis import (
+    certify_jerk_cost,
+    certify_limits,
+    check_limits,
+    compute_metrics,
+)
 from trajectory_optimizer.model import (
     BoundaryConditions,
     BoundaryKind,
@@ -37,6 +43,7 @@
     "BoundaryConditions",
     "BoundaryKind",
     "CertifiedTrajectory",
+    "IntegrableTrajectory",
     "JointLimits",
     "LimitSet",
     "PeakDerivatives",
@@ -49,6 +56,7 @@
     "ViaPoints",
     "Waypoint",
     "__version__",
+    "certify_jerk_cost",
     "certify_limits",
     "check_limits",
     "compute_metrics",

src/trajectory_optimizer/algorithm/__init__.py

@@ -20,6 +20,7 @@
 from trajectory_optimizer.algorithm.protocol import (
     DERIVATIVE_NAMES,
     CertifiedTrajectory,
+    IntegrableTrajectory,
     PeakDerivatives,
     Trajectory,
     TrajectoryState,
@@ -50,6 +51,7 @@
 __all__ = [
     "DERIVATIVE_NAMES",
     "CertifiedTrajectory",
+    "IntegrableTrajectory",
     "LimitActivation",
     "PeakDerivatives",
     "PiecewisePolynomialTrajectory",

src/trajectory_optimizer/algorithm/piecewise.py

@@ -230,6 +230,35 @@ def peak_derivatives(self) -> PeakDerivatives:
                         times[row, joint] = origin + float(candidates[best])
         return PeakDerivatives(values=values, times=times)
 
+    def jerk_cost(self) -> float:
+        """Integrated squared jerk summed over joints, solved in closed form.
+
+        On one segment the jerk of one joint is a polynomial in the local time,
+        so its square is a polynomial too and the integral over that segment is
+        the quadratic form ``c @ M @ c`` in the jerk coefficients ``c``, with
+
+        ``M[i, k] = span ** (i + k + 1) / (i + k + 1)``.
+
+        Summing that over every segment and every joint gives the objective the
+        minimum-jerk generator minimises, with no quadrature error: the value is
+        exact rather than converging as a sample count grows. That matters
+        because the integrand of a jerk limited profile is piecewise constant
+        and discontinuous, and Simpson's rule then converges only as the
+        reciprocal of the sample count.
+
+        A profile whose acceleration steps at a breakpoint, the trapezoidal one
+        for instance, carries an impulse of jerk there that no integral of the
+        segment polynomials can see, so the cost reported for it is zero. That
+        is the same one-sided reading of a piecewise definition that
+        :meth:`peak_derivatives` gives, and it is not a bound either.
+        """
+        table = self._derivatives[3]
+        orders = np.arange(table.shape[2], dtype=np.float64)
+        exponents = orders[:, None] + orders[None, :] + 1.0
+        spans = np.diff(self._breakpoints)[:, None, None]
+        weights = spans**exponents / exponents
+        return float(np.einsum("sji,sik,sjk->", table, weights, table))
+
     def time_scaled(self, factor: float) -> PiecewisePolynomialTrajectory:
         """Return the same geometric path traversed ``factor`` times more slowly.
 

src/trajectory_optimizer/algorithm/protocol.py

@@ -16,6 +16,7 @@
 __all__ = [
     "DERIVATIVE_NAMES",
     "CertifiedTrajectory",
+    "IntegrableTrajectory",
     "PeakDerivatives",
     "Trajectory",
     "TrajectoryState",
@@ -111,6 +112,23 @@ def peak_derivatives(self) -> PeakDerivatives:
         ...
 
 
+@runtime_checkable
+class IntegrableTrajectory(Protocol):
+    """A trajectory that integrates its own squared jerk in closed form.
+
+    A quadrature over a sampled trace converges to the jerk cost only as fast as
+    the integrand is smooth, and the integrand of a jerk limited profile is
+    discontinuous. A trajectory that implements this protocol integrates its own
+    coefficient table instead, which carries no quadrature error at all.
+    Wrappers raise :class:`TypeError` when the trajectory they wrap cannot do
+    the same.
+    """
+
+    def jerk_cost(self) -> float:
+        """Exact integral of the squared jerk, summed over joints."""
+        ...
+
+
 @dataclass(frozen=True, slots=True, eq=False)
 class TrajectoryState:
     """Position and its first three derivatives at a single instant."""

src/trajectory_optimizer/algorithm/scaling.py

@@ -20,6 +20,7 @@
 from trajectory_optimizer.algorithm.protocol import (
     DERIVATIVE_NAMES,
     CertifiedTrajectory,
+    IntegrableTrajectory,
     PeakDerivatives,
     Trajectory,
 )
@@ -111,6 +112,20 @@ def peak_derivatives(self) -> PeakDerivatives:
             times=peaks.times * self._scale,
         )
 
+    def jerk_cost(self) -> float:
+        """Exact jerk cost of the base, divided by the fifth power of the scale.
+
+        Stretching time divides the jerk by ``scale ** 3`` and multiplies the
+        interval it is integrated over by ``scale``, so the integral of the
+        squared jerk is divided by ``scale ** 5``. Raises :class:`TypeError`
+        when the base cannot integrate its own squared jerk in closed form.
+        """
+        if not isinstance(self._base, IntegrableTrajectory):
+            raise TypeError(
+                f"{type(self._base).__name__} cannot integrate its squared jerk in closed form"
+            )
+        return self._base.jerk_cost() / self._scale**5
+
 
 @dataclass(frozen=True, slots=True)
 class LimitActivation:

src/trajectory_optimizer/algorithm/synchronisation.py

@@ -17,6 +17,7 @@
 from trajectory_optimizer.algorithm.piecewise import constant_trajectory
 from trajectory_optimizer.algorithm.protocol import (
     CertifiedTrajectory,
+    IntegrableTrajectory,
     PeakDerivatives,
     Trajectory,
 )
@@ -98,6 +99,24 @@ def peak_derivatives(self) -> PeakDerivatives:
             times=np.concatenate([block.times for block in blocks], axis=1),
         )
 
+    def jerk_cost(self) -> float:
+        """Exact jerk cost of every component, added together.
+
+        The cost is a sum over joints and each component owns a disjoint block
+        of them, so adding the component costs is exact. Raises
+        :class:`TypeError` when any component cannot integrate its own squared
+        jerk in closed form.
+        """
+        total = 0.0
+        for index, component in enumerate(self._components):
+            if not isinstance(component, IntegrableTrajectory):
+                raise TypeError(
+                    f"component {index} of type {type(component).__name__} cannot integrate "
+                    "its squared jerk in closed form"
+                )
+            total += component.jerk_cost()
+        return total
+
 
 def synchronise(components: Sequence[Trajectory]) -> SynchronisedTrajectory:
     """Stretch every component onto the duration of the slowest one.

src/trajectory_optimizer/analysis/__init__.py

@@ -12,7 +12,11 @@
     certify_limits,
     check_limits,
 )
-from trajectory_optimizer.analysis.metrics import TrajectoryMetrics, compute_metrics
+from trajectory_optimizer.analysis.metrics import (
+    TrajectoryMetrics,
+    certify_jerk_cost,
+    compute_metrics,
+)
 from trajectory_optimizer.analysis.report import (
     ProfileSummary,
     format_comparison_table,
@@ -29,6 +33,7 @@
     "LimitViolation",
     "ProfileSummary",
     "TrajectoryMetrics",
+    "certify_jerk_cost",
     "certify_limits",
     "check_limits",
     "comparison_figure",