Skip to content

Commit e205500

Browse files
committed
feat(motion_profiles): add trapezoidal (LSPB) velocity profile
Deploy roadmap item 56: minimum-time trapezoid / triangle plan from (vMax, aMax), plus PlanWithDuration, which keeps aMax and lowers the cruise speed so several axes can finish together. The cruise speed uses the cancellation-free conjugate form of the smaller root, and infeasible durations return nullopt. At and after tf the profile holds the goal at rest; before t = 0 it holds the start. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Q53m7MZDmHwocmKbhhd2o5
1 parent 80f2861 commit e205500

13 files changed

Lines changed: 355 additions & 226 deletions

File tree

‎README.md‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,7 @@ Refer to the documentation to quickly integrate and utilize the library's signal
2424
| [Optimization](doc/optimization/README.md) | Gradient Descent, SGD / Adam Step Optimisers |
2525
| [Regularization](doc/regularization/README.md) | L1 (Lasso), L2 (Ridge) |
2626
| [Math](doc/math/README.md) | CORDIC, Quaternion, MatrixNorms, Step Response Metrics, MatrixExponential, Angle Utilities, Tridiagonal (Thomas) Solver |
27-
| [Motion Profiles](doc/motion_profiles/README.md) | Polynomial (Cubic/Quintic) Point-to-Point Profile |
27+
| [Motion Profiles](doc/motion_profiles/README.md) | Polynomial (Cubic/Quintic) Point-to-Point Profile, Trapezoidal (LSPB) Profile |
2828
| [Solvers](doc/solvers/README.md) | Gaussian Elimination, Levinson-Durbin, Durand-Kerner, Cholesky, DARE, Runge-Kutta ODE Integrators (RK4 + Dormand-Prince), Spectral Radius & Discrete Stability Margin, QR Decomposition (Householder / Givens), LU Decomposition with Partial Pivoting, Singular Value Decomposition (Golub-Kahan), Symplectic Integrators (Semi-Implicit Euler, Velocity Verlet) |
2929
| [Nonlinear Control](doc/nonlinear_control/README.md) | Feedback Linearization, Backstepping Control, Model Reference Adaptive Control (MRAC) |
3030
| [Robust Control](doc/robust_control/README.md) | Active Disturbance Rejection Control (ADRC + ESO), Sliding Mode Control (SMC), Disturbance Observer (DOB), H∞ State-Feedback Control |

‎ROADMAP.md‎

Lines changed: 0 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,6 @@ Difficulty legend:
2626

2727
| # | Component | Target module | Difficulty |
2828
|----|------------------------------------------------------|-------------------------|------------|
29-
| 56 | Trapezoidal (LSPB) profile + duration sync | `motion_profiles` (new) | ★☆☆☆☆ |
3029
| 53 | SE(3) transform, twists, wrenches, adjoint, exp/log | `math` | ★★☆☆☆ |
3130
| 60 | S-curve (jerk-limited) profile | `motion_profiles` (new) | ★★☆☆☆ |
3231
| 61 | Cubic spline interpolation (vector-valued) | `math` | ★★☆☆☆ |
@@ -44,17 +43,6 @@ unit-test reference [`TESTING.md`](TESTING.md) depends on but which
4443
the library does not yet expose. Detailed below under
4544
[Evaluation & metrics primitives](#evaluation--metrics-primitives).
4645

47-
## Tier 1 — Trivial primitives ★☆☆☆☆
48-
49-
### 56. Trapezoidal (LSPB) profile
50-
- **What:** Minimum-time velocity/acceleration-limited rest-to-rest move, triangular for short moves, plus
51-
a fixed-duration plan for multi-axis synchronization.
52-
- **Embedded value:** The standard servo-axis move (robot joints, CNC, conveyors, motor position loops).
53-
- **Algorithm / paper:** L. Biagiotti, C. Melchiorri, *Trajectory Planning for Automatic Machines and
54-
Robots* (2008), Ch. 3.
55-
- **Reuses:** Item 55 types.
56-
- **Spec:** [roadmap/motion_profiles/TrapezoidalProfile](roadmap/motion_profiles/TrapezoidalProfile/implementation.md).
57-
5846
## Tier 2 — Easy ★★☆☆☆
5947

6048
### 53. SE(3) transform, twists and wrenches

‎doc/motion_profiles/README.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,3 +7,4 @@ Time laws for point-to-point motion of a single axis: position, velocity, accele
77
| Algorithm | Description |
88
|-------------------------------------------------|---------------------------------------------------------------------------------------------|
99
| [Polynomial Profile](PolynomialProfile.md) | Cubic (position + velocity) or quintic (+ acceleration) boundary-matched time law over a fixed duration |
10+
| [Trapezoidal Profile](TrapezoidalProfile.md) | Minimum-time trapezoid / triangle velocity profile (LSPB) and a fixed-duration plan for multi-axis synchronization |
Lines changed: 111 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,111 @@
1+
# Trapezoidal Profile (Linear Segment with Parabolic Blends)
2+
3+
## Overview & Motivation
4+
5+
A point-to-point move has to reach its goal as fast as the actuator allows without saturating the drive. The trapezoidal profile — *linear segment with parabolic blends* (LSPB) — does this with the simplest rule that respects both a velocity ceiling and an acceleration ceiling: accelerate at the maximum rate, coast at the maximum speed, decelerate at the maximum rate. Velocity plotted against time is a trapezoid; a move too short to reach the ceiling becomes a triangle.
6+
7+
From the limits and the distance alone, the plan is the **minimum-time** one. A second mode fixes the duration instead, which multi-axis synchronization needs: all axes must finish together, so the faster ones are slowed rather than the slowest sped up.
8+
9+
## Mathematical Theory
10+
11+
### Problem Statement
12+
13+
Rest-to-rest move of length $d = |q_f - q_0|$ in direction $\sigma = \pm 1$ (positive when $q_f \ge q_0$), with limits $v_{max} > 0$ and $a_{max} > 0$. The profile is planned for the magnitude and mirrored by $\sigma$.
14+
15+
### Minimum-Time Plan
16+
17+
A ramp between rest and speed $v$ takes $v/a_{max}$ and covers $v^2/(2a_{max})$, so the two ramps together consume the **blend distance** $v^2/a_{max}$. Comparing $d$ with the blend distance at $v_{max}$ selects the shape:
18+
19+
$$
20+
\begin{cases}
21+
d \ge v_{max}^2/a_{max}: & v = v_{max}, \quad t_a = v/a_{max}, \quad t_c = (d - v^2/a_{max})/v \quad\text{(trapezoid)} \\
22+
d < v_{max}^2/a_{max}: & v = \sqrt{d\,a_{max}}, \quad t_a = v/a_{max}, \quad t_c = 0 \quad\text{(triangle)}
23+
\end{cases}
24+
$$
25+
26+
The total time is $T = 2t_a + t_c$. The cases agree at $d = v_{max}^2/a_{max}$, so $T$ is continuous in $d$.
27+
28+
### Sampling
29+
30+
With $v$, $t_a$, $t_c$, $T$ precomputed, the unsigned state is:
31+
32+
| Phase | Interval | Acceleration | Velocity | Position |
33+
|-------|----------|--------------|----------|----------|
34+
| Parabolic ramp-up | $0 \le t < t_a$ | $a_{max}$ | $a_{max}\,t$ | $\tfrac12 a_{max}\,t^2$ |
35+
| Linear cruise | $t_a \le t < t_a + t_c$ | $0$ | $v$ | $v\,(t - t_a/2)$ |
36+
| Parabolic ramp-down | $t_a + t_c \le t < T$ | $-a_{max}$ | $a_{max}\,(T-t)$ | $d - \tfrac12 a_{max}\,(T-t)^2$ |
37+
38+
Velocity, acceleration and the position offset are multiplied by $\sigma$. At and beyond $t = T$ the profile holds the goal at rest (velocity and acceleration zero); times before $0$ are treated as $t = 0$, holding the start with zero velocity. The reported jerk is zero in every phase.
39+
40+
### Fixed-Duration Plan
41+
42+
The acceleration $a_{max}$ is kept and the cruise speed is lowered until the move takes exactly $T$. The ramps cover $v\,t_a$ in total and the cruise covers $v\,(T - 2t_a)$, so with $t_a = v/a_{max}$
43+
44+
$$
45+
d = v\left(T - \frac{v}{a_{max}}\right) \;\Longleftrightarrow\; v^2 - a_{max}T\,v + a_{max}\,d = 0 .
46+
$$
47+
48+
The smaller root keeps $t_c \ge 0$. Multiplying it by its conjugate gives a form with no subtraction of nearly equal numbers:
49+
50+
$$
51+
v = \frac{a_{max}T - \sqrt{a_{max}^2T^2 - 4a_{max}d}}{2} = \frac{2\,a_{max}\,d}{a_{max}T + \sqrt{a_{max}^2T^2 - 4a_{max}d}} .
52+
$$
53+
54+
A real root exists iff $a_{max}\,T^2 \ge 4d$. Equality gives a double root $v = a_{max}T/2$ and a triangle: $T$ is then the shortest duration achievable at this acceleration. Otherwise $t_a = v/a_{max}$ and $t_c = T - 2t_a$.
55+
56+
### Multi-Axis Synchronization
57+
58+
Let $T_{sync}$ be the longest axis time and $\lambda = T_{sync}/T_{axis} \ge 1$. Planning an axis again with $(v_{max}/\lambda,\; a_{max}/\lambda^2)$ leaves the blend distance $v_{max}^2/a_{max}$ unchanged, so the same shape is selected, while $t_a$ and $t_c$ both scale by $\lambda$: the duration is exactly $\lambda\,T_{axis}$. Velocity is divided by $\lambda$ and acceleration by $\lambda^2$.
59+
60+
The fixed-duration plan is the alternative. It keeps $a_{max}$ and only lowers the cruise speed: the smaller root decreases as $T$ grows, so for $T \ge T_{min}$ the cruise speed never exceeds $v_{max}$.
61+
62+
## Complexity Analysis
63+
64+
| Operation | Time | Space | Notes |
65+
|-----------|------|-------|-------|
66+
| Minimum-time plan | $O(1)$ | $O(1)$ | One comparison, one square root (triangle only) |
67+
| Fixed-duration plan | $O(1)$ | $O(1)$ | One square root, one division |
68+
| Sample | $O(1)$ | $O(1)$ | One phase branch, a few multiply-adds; no state change |
69+
70+
Memory is a fixed handful of scalars: start, direction, distance, acceleration, peak speed and three phase times.
71+
72+
## Step-by-Step Walkthrough
73+
74+
**Trapezoid.** $d = 5$, $v_{max} = 1$, $a_{max} = 2$. The blend distance is $1/2 = 0.5 \le 5$, so a cruise exists: $t_a = 0.5$, $t_c = (5 - 0.5)/1 = 4.5$, $T = 2(0.5) + 4.5 = 5.5$.
75+
76+
| $t$ | Phase | Position | Velocity | Acceleration |
77+
|-----|-------|----------|----------|--------------|
78+
| $0.25$ | ramp-up | $\tfrac12\cdot 2\cdot 0.25^2 = 0.0625$ | $0.5$ | $2$ |
79+
| $3$ | cruise | $1\cdot(3 - 0.25) = 2.75$ | $1$ | $0$ |
80+
| $5.25$ | ramp-down | $5 - \tfrac12\cdot 2\cdot 0.25^2 = 4.9375$ | $0.5$ | $-2$ |
81+
| $\ge 5.5$ | held | $5$ | $0$ | $0$ |
82+
83+
The ramp-up ends at position $0.25$, and the cruise formula gives the same value at $t = 0.5$. The triangle threshold for these limits is $d = 0.5$; for $d = 0.25$, $v = \sqrt{0.5} \approx 0.7071$ and $T \approx 0.7071$, with no cruise.
84+
85+
**Fixed duration.** $d = 4$, $a_{max} = 1$, $T = 5$: discriminant $25 - 16 = 9$, $v = 2\cdot 4/(5 + 3) = 1$, $t_a = 1$, $t_c = 3$; check $1\cdot(5 - 1) = 4$. For $d = 5$, $a_{max} = 2$, $T = 6$: discriminant $144 - 40 = 104$, $v = 20/(12 + \sqrt{104}) \approx 0.900980$.
86+
87+
## Pitfalls & Edge Cases
88+
89+
- **Non-positive limits.** The minimum-time plan requires $v_{max} > 0$ and $a_{max} > 0$; a violation is an assertion failure, not a recoverable status. The fixed-duration plan returns no result for $a_{max} \le 0$ or $T \le 0$.
90+
- **Infeasible duration.** If $a_{max}T^2 < 4d$ the plan is rejected, never clamped; the shortest feasible duration is $2\sqrt{d/a_{max}}$. At the boundary, single-precision rounding can push the discriminant just below zero, so plan the slowest axis with the minimum-time plan.
91+
- **Cancellation.** The textbook root $(a_{max}T - \sqrt{\cdot})/2$ loses its significant digits in single precision when $4a_{max}d \ll a_{max}^2T^2$; the conjugate form does not. The computed cruise time is clamped at zero.
92+
- **Zero distance.** The minimum-time plan has $v = 0$ and $T = 0$: a zero-length profile that is at rest at the goal for every $t \ge 0$.
93+
- **Discontinuous acceleration.** Position and velocity are continuous, but acceleration jumps between $0$ and $\pm a_{max}$ at the joins. The jerk impulses are not represented in the reported jerk and can excite structural modes; use the S-curve profile in that case.
94+
95+
## Variants & Generalizations
96+
97+
The **S-curve** profile adds a jerk limit (seven segments) for continuous acceleration at the cost of a longer move. **Polynomial** profiles (cubic, quintic) are fixed-duration and smooth by construction but enforce no velocity or acceleration ceiling. Unequal acceleration and deceleration limits use the same construction with unequal ramp times.
98+
99+
## Applications
100+
101+
- **Point-to-point axis moves** — servo and stepper drives, gantries, pick-and-place, gimbal slews.
102+
- **Multi-axis coordination** — time-scaled or fixed-duration plans make every axis arrive together.
103+
104+
## Connections to Other Algorithms
105+
106+
The S-curve profile is the jerk-limited upgrade under the same velocity and acceleration limits. The polynomial profile is fixed-time and smooth but offers no limit guarantees. The saturation rate limiter slews a changing setpoint online, whereas the trapezoid plans a whole move offline.
107+
108+
## References & Further Reading
109+
110+
- Biagiotti, L. and Melchiorri, C., "Trajectory Planning for Automatic Machines and Robots," Springer, 2008, §3.2.
111+
- Siciliano, B., Sciavicco, L., Villani, L., and Oriolo, G., "Robotics: Modelling, Planning and Control," Springer, 2009, §4.2.

‎numerical/motion_profiles/CMakeLists.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -13,10 +13,12 @@ target_link_libraries(numerical.motion_profiles ${NUMERICAL_VISIBILITY}
1313
target_sources(numerical.motion_profiles PRIVATE
1414
MotionProfileTypes.hpp
1515
PolynomialProfile.hpp
16+
TrapezoidalProfile.hpp
1617
)
1718

1819
numerical_add_coverage_sources(numerical.motion_profiles
1920
PolynomialProfile.cpp
21+
TrapezoidalProfile.cpp
2022
)
2123

2224
add_subdirectory(test)
Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,6 @@
1+
#include "numerical/motion_profiles/TrapezoidalProfile.hpp"
2+
3+
namespace motion_profiles
4+
{
5+
template class TrapezoidalProfile<float>;
6+
}
Lines changed: 114 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,114 @@
1+
#pragma once
2+
3+
#include "infra/util/ReallyAssert.hpp"
4+
#include "numerical/math/CompilerOptimizations.hpp"
5+
#include "numerical/math/Math.hpp"
6+
#include "numerical/motion_profiles/MotionProfileTypes.hpp"
7+
#include <algorithm>
8+
#include <optional>
9+
#include <type_traits>
10+
11+
namespace motion_profiles
12+
{
13+
template<typename T>
14+
class TrapezoidalProfile
15+
{
16+
static_assert(std::is_floating_point_v<T>, "TrapezoidalProfile supports floating-point types only");
17+
18+
public:
19+
TrapezoidalProfile(T q0, T qf, const MotionLimits<T>& limits);
20+
static std::optional<TrapezoidalProfile> PlanWithDuration(T q0, T qf, T aMax, T duration);
21+
22+
OPTIMIZE_FOR_SPEED TrajectoryState<T> Sample(T t) const;
23+
T Duration() const;
24+
bool IsTriangular() const;
25+
26+
private:
27+
TrapezoidalProfile(T q0, T qf, T aMax);
28+
29+
T q0;
30+
T direction;
31+
T distance;
32+
T aMax;
33+
T vPeak{};
34+
T tAccel{};
35+
T tCruise{};
36+
T tf{};
37+
};
38+
39+
template<typename T>
40+
TrapezoidalProfile<T>::TrapezoidalProfile(T q0, T qf, T aMax)
41+
: q0{ q0 }
42+
, direction{ qf >= q0 ? T{ 1 } : T{ -1 } }
43+
, distance{ math::Abs(qf - q0) }
44+
, aMax{ aMax }
45+
{}
46+
47+
template<typename T>
48+
TrapezoidalProfile<T>::TrapezoidalProfile(T q0, T qf, const MotionLimits<T>& limits)
49+
: TrapezoidalProfile{ q0, qf, limits.aMax }
50+
{
51+
really_assert(limits.vMax > T{});
52+
really_assert(limits.aMax > T{});
53+
54+
const T blendDistance{ limits.vMax * limits.vMax / limits.aMax };
55+
if (distance >= blendDistance)
56+
{
57+
vPeak = limits.vMax;
58+
tCruise = (distance - blendDistance) / vPeak;
59+
}
60+
else
61+
vPeak = math::Sqrt(distance * aMax);
62+
63+
tAccel = vPeak / aMax;
64+
tf = T{ 2 } * tAccel + tCruise;
65+
}
66+
67+
template<typename T>
68+
std::optional<TrapezoidalProfile<T>> TrapezoidalProfile<T>::PlanWithDuration(T q0, T qf, T aMax, T duration)
69+
{
70+
TrapezoidalProfile profile{ q0, qf, aMax };
71+
const T discriminant{ aMax * aMax * duration * duration - T{ 4 } * aMax * profile.distance };
72+
if (aMax <= T{} || duration <= T{} || discriminant < T{})
73+
return std::nullopt;
74+
75+
profile.vPeak = T{ 2 } * aMax * profile.distance / (aMax * duration + math::Sqrt(discriminant));
76+
profile.tAccel = profile.vPeak / aMax;
77+
profile.tCruise = std::max(duration - T{ 2 } * profile.tAccel, T{});
78+
profile.tf = duration;
79+
return profile;
80+
}
81+
82+
template<typename T>
83+
OPTIMIZE_FOR_SPEED TrajectoryState<T> TrapezoidalProfile<T>::Sample(T t) const
84+
{
85+
t = std::max(t, T{});
86+
if (t >= tf)
87+
return { q0 + direction * distance, T{}, T{}, T{} };
88+
89+
if (t < tAccel)
90+
return { q0 + direction * T{ 0.5 } * aMax * t * t, direction * aMax * t, direction * aMax, T{} };
91+
92+
if (t < tAccel + tCruise)
93+
return { q0 + direction * vPeak * (t - T{ 0.5 } * tAccel), direction * vPeak, T{}, T{} };
94+
95+
const T remaining{ tf - t };
96+
return { q0 + direction * (distance - T{ 0.5 } * aMax * remaining * remaining), direction * aMax * remaining, -direction * aMax, T{} };
97+
}
98+
99+
template<typename T>
100+
T TrapezoidalProfile<T>::Duration() const
101+
{
102+
return tf;
103+
}
104+
105+
template<typename T>
106+
bool TrapezoidalProfile<T>::IsTriangular() const
107+
{
108+
return tCruise <= T{};
109+
}
110+
111+
#ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD
112+
extern template class TrapezoidalProfile<float>;
113+
#endif
114+
}

‎numerical/motion_profiles/test/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,4 +11,5 @@ numerical_link_qemu_runtime(numerical.motion_profiles_test)
1111

1212
target_sources(numerical.motion_profiles_test PRIVATE
1313
TestPolynomialProfile.cpp
14+
TestTrapezoidalProfile.cpp
1415
)

0 commit comments

Comments
 (0)