Skip to content

Latest commit

 

History

History
1239 lines (870 loc) · 42.2 KB

File metadata and controls

1239 lines (870 loc) · 42.2 KB

Algorithms

Technical reference for every algorithm used in fullstack-manip. Each section documents the theoretical background, the mathematical formulation, and the code implementation.


Table of Contents


Part I — Kinematics & Planning

This part covers algorithms that determine where the robot should be (inverse kinematics), how it gets there (trajectory planning), and where it can reach (workspace sampling).


1.1 Differential IK (Mink/DAQP)

Theory

Differential IK solves the first-order linearised kinematics:

$$\dot{q} = J^+ \dot{x}$$

where $J \in \mathbb{R}^{6 \times n}$ is the geometric Jacobian of the end-effector frame. The Mink library formulates this as a Quadratic Programme (QP) at each integration step:

$$\min_{\dot{q}} \quad \tfrac{1}{2}|\dot{q}|^2 \quad \text{s.t.} \quad J\dot{q} = \dot{x}^*, \quad \dot{q}^{lo} \leq \dot{q} \leq \dot{q}^{hi}$$

The position/orientation task error is:

$$e = \begin{bmatrix} p^* - p \ \log(R^* R^{-1})^\vee \end{bmatrix} \in \mathbb{R}^6$$

The QP is solved by the DAQP active-set solver.

Each joint-velocity solution is integrated forward with Euler:

$$q_{k+1} = q_k + \dot{q}_k \cdot \Delta t$$

Convergence thresholds and iteration cap are configurable constructor kwargs (defaults: pos_threshold=1e-4, ori_threshold=1e-4, max_iters=100). Custom IK backends can be injected via ik_fn — all Mink/MuJoCo logic is bypassed when provided.

Parameters

Name Value Description
SOLVER "daqp" QP backend
pos_threshold 1e-4 m Position convergence (configurable)
ori_threshold 1e-4 rad Orientation convergence (configurable)
max_iters 100 Max integration steps (configurable)
dt 0.01 s Integration timestep
position_cost 1.0 Task weight
lm_damping 1e-6 Levenberg-Marquardt regularisation

Implementation

fullstack_manip/kinematics/ik_solver.py
  IKSolver
    solve_ik_for_qpos(frame_name, frame_type, target_pos, target_orient) → q[6]

Key steps:

# 1. Build SE(3) target
target = mink.SE3.from_rotation_and_translation(SO3(target_orient), target_pos)

# 2. Create frame task
task = mink.FrameTask(frame_name, frame_type,
                      position_cost=1.0, orientation_cost=0.0, lm_damping=1e-6)
task.set_target(target)

# 3. Integrate until convergence
for i in range(self._max_iters):
    vel = mink.solve_ik(config, [task], dt, solver="daqp", limits=limits)
    config.integrate_inplace(vel, dt)
    e = task.compute_error(config)
    if norm(e[:3]) <= self._pos_threshold and norm(e[3:]) <= self._ori_threshold:
        break

return config.q[:6]

1.2 Joint-Space Cubic Spline Planning

Theory

Given start and goal Cartesian poses $(x_0, x_f)$, the planner:

  1. IK → maps each Cartesian pose to a joint configuration:

$$q_0 = \text{IK}(x_0), \qquad q_f = \text{IK}(x_f)$$

  1. Duration heuristic — minimum time to respect maximum joint velocity:

$$T = \max!\left(\frac{|q_f - q_0|_\infty}{v_{max}},\ T_{desired}\right)$$

  1. Cubic spline$C^2$-continuous trajectory interpolating boundary conditions:

$$q(t) = \text{CubicSpline}!\bigl([0,, T],\ [q_0,, q_f]\bigr)$$

sampled at uniform timesteps $t_k = k,\Delta t$, producing a $(T/\Delta t,, n_q)$ array.

The cubic polynomial per segment satisfies:

$$q(0) = q_0,\quad q(T) = q_f,\quad \dot{q}(0) = v_0,\quad \dot{q}(T) = v_f$$

where $v_0, v_f$ default to zero (rest-to-rest motion).

Parameters

Name Default Description
max_velocity 0.5 rad/s Per-joint velocity bound
max_acceleration 1.0 rad/s² Per-joint acceleration bound
dt 0.01 s Trajectory sample rate

Implementation

fullstack_manip/planning/motion_planner.py
  MotionPlanner
    plan_trajectory(start_pos, end_pos, ...) → (trajectory[T,nq], t[T])
    solve_ik_for_qpos(frame_name, ...) → q[6]

The planner delegates IK to IKSolver, trajectory generation to scipy.interpolate.CubicSpline, and collision checking to the injected CollisionChecker instance.


1.3 Reachability Workspace Sampling

Theory

The reachable workspace $\mathcal{W}$ is estimated by Monte Carlo grid sampling: a 3-D grid of candidate end-effector positions is generated around a centre point, IK is solved for each, and the positions where IK succeeds and the FK error is below tolerance are retained.

Algorithm:

  1. Grid: axis = linspace(-R, R, G) along each axis, full meshgrid → $G^3$ candidates.
  2. Random subsample to $M$ using np.random.default_rng(seed) for reproducibility.
  3. For each candidate $\mathbf{p}^* = \mathbf{c} + \mathbf{offset}$:
    • Solve $q = \text{IK}(\mathbf{p}^*)$ from current seed configuration $q_{seed}$
    • Apply $q$, run FK, get $\mathbf{p}_{fk}$
    • Accept if $|\mathbf{p}{fk} - \mathbf{p}^*| \leq \delta{tol}$
    • Update $q_{seed} \leftarrow q$ (warm-start for nearby points)
  4. Deduplicate by rounding to $\delta_{tol}$ and np.unique.

Parameters

Name Default Description
grid_size $G$ 20 Points per axis ($G^3$ total candidates)
grid_range $R$ 0.4 m Half-width of the sampling volume
max_points $M$ 4000 Subsample cap for IK budget
tolerance $\delta_{tol}$ 5e-3 m FK acceptance threshold
seed 0 RNG seed for reproducibility

Implementation

fullstack_manip/core/robot.py
  Robot
    compute_workspace(ik_solver, grid_size=20, grid_range=0.4,
                      max_points=4000, tolerance=5e-3, center, seed=0)
                      → np.ndarray[(N,3)]
    visualize_workspace(points)   # renders mjGEOM_SPHERE via user_scn

Part II — Control

This part covers algorithms that execute motion: low-level joint regulation, Cartesian-space interaction controllers (impedance, admittance, force), model-based gravity/friction cancellation, and the pick-and-place task sequencer.


2.1 PID Joint Tracking

Theory

Each joint is regulated by an independent discrete-time PID controller:

$$u_k = K_p, e_k + K_i\sum_{j=0}^{k} e_j,\Delta t + K_d,\frac{e_k - e_{k-1}}{\Delta t}$$

with an additional passive damping term ($K_{damp} = 0.1$) subtracted from the derivative component to reduce high-frequency chatter:

$$u_k = K_p, e_k + K_i \int e, dt + (K_d - 0.1),\dot{e}$$

Optional clamping: integral is clipped to $\pm I_{max}$; output is clipped to $\pm u_{max}$.

The MotionExecutor wraps the PID loop with a contact-stop condition: trajectory execution halts when CollisionChecker.detect_contact(ee, object_geom) returns True.

Default Gains

Gain Default
$K_p$ 2.0 (per joint, uniform)
$K_i$ 0.1
$K_d$ 0.5

Implementation

fullstack_manip/control/motion_executor.py
  MotionExecutor
    move_to_pose(target_pos, target_orient, duration, ...)

fullstack_manip/control/low_level/pid_controller.py
  PIDController
    compute(target, current, timestamp=None) → u[n_q]
    reset()
    set_gains(kp, ki, kd)

2.2 Cartesian Impedance Control

Theory

The robot is modelled as a virtual mass-spring-damper in Cartesian space. The impedance relationship between end-effector displacement and applied force is:

$$\mathbf{F}_{cmd} = \mathbf{M},\ddot{\mathbf{x}} + \mathbf{B},\dot{\mathbf{x}} + \mathbf{K},\Delta\mathbf{x}$$

where $\Delta\mathbf{x} = \mathbf{x}^* - \mathbf{x} \in \mathbb{R}^6$ is the 6-DOF Cartesian error (3 position + 3 orientation).

The current implementation uses the proportional (stiffness) term in the motion loop:

$$\mathbf{F}_{cmd} = \mathbf{K}_{3\times3},\Delta\mathbf{p}$$ $$\delta\mathbf{p} = \mathbf{F}_{cmd},\Delta t, \qquad \mathbf{p}_{new} = \mathbf{p} + \delta\mathbf{p}$$

Design rationale: all three matrices $\mathbf{M}, \mathbf{B}, \mathbf{K} \in \mathbb{R}^{6\times6}$ must be symmetric positive-definite (SPD). The constructor enforces this by allowing scalar ($\rightarrow k\mathbf{I}_6$), 6-vector ($\rightarrow \text{diag}$), or full $6\times6$ input.

Critical damping condition: $\mathbf{B} = 2\sqrt{\mathbf{K}\mathbf{M}}$ (component-wise for diagonal matrices).

Parameters

Matrix Default Unit
$\mathbf{K}$ $500,\mathbf{I}_6$ N/m, Nm/rad
$\mathbf{B}$ $50,\mathbf{I}_6$ Ns/m, Nms/rad
$\mathbf{M}$ $1.0,\mathbf{I}_6$ kg, kg·m²

Implementation

fullstack_manip/control/high_level/impedance_controller.py
  ImpedanceController(BaseController)
    execute(desired_pos, desired_orient, duration=5.0)
    set_stiffness(value, axis=None)   # axis ∈ {x,y,z,rx,ry,rz} → index 0–5
    set_damping(value, axis=None)
    compliant_approach(target_pos, low_stiffness=100.0)
    absorb_impact(duration=2.0, damping_ratio=2.0)
    get_impedance_parameters() → {'stiffness': K, 'damping': B, 'inertia': M}

2.3 Cartesian Admittance Control

Theory

Admittance control is the dual of impedance: external forces are mapped to commanded Cartesian velocities. The admittance model is:

$$\mathbf{M},\ddot{\mathbf{x}} + \mathbf{B},\dot{\mathbf{x}} + \mathbf{K},(\mathbf{x} - \mathbf{x}_{ref}) = \mathbf{F}_{ext}$$

Rearranging for acceleration and integrating with forward Euler:

$$\ddot{\mathbf{x}}_k = \mathbf{M}^{-1}!\left(\mathbf{F}_{ext} - \mathbf{B},\dot{\mathbf{x}}_k - \mathbf{K},(\mathbf{x}_k - \mathbf{x}_{ref})\right)$$ $$\dot{\mathbf{x}}_{k+1} = \dot{\mathbf{x}}_k + \ddot{\mathbf{x}}_k,\Delta t$$ $$\mathbf{x}_{k+1} = \mathbf{x}_k + \dot{\mathbf{x}}_k,\Delta t$$

The linear system $\mathbf{M},\ddot{\mathbf{x}} = \cdot$ is solved via np.linalg.solve (LU factorisation of $\mathbf{M}$).

Stability condition: $\mathbf{M}, \mathbf{B} \succ 0$. Pure admittance is obtained by setting $\mathbf{K} = \mathbf{0}$.

Velocity decay: when $|\mathbf{F}{ext}[:3]| < F{threshold}$, velocity is multiplied by 0.95 per timestep (exponential decay to rest).

External force $\mathbf{F}_{ext}$ is read from data.contact[i].frame[:3] and accumulated into a 6-vector (3 force + 3 torque).

Parameters

Matrix Default Description
$\mathbf{M}$ $5.0,\mathbf{I}_6$ Virtual inertia
$\mathbf{B}$ $20.0,\mathbf{I}_6$ Virtual damping
$\mathbf{K}$ $\mathbf{0}$ Stiffness (zero = pure admittance)

Implementation

fullstack_manip/control/high_level/admittance_controller.py
  AdmittanceController(BaseController)
    execute(duration=10.0, force_threshold=1.0)
    follow_force(force_source, duration, velocity_limit=0.1)
    backdrive_motion(duration=5.0, sensitivity=0.5) → [(pos, orient)]
    kinesthetic_teaching(duration=10.0) → [q]
    _compute_acceleration(F_ext) → ẍ[6]
    _measure_external_force() → F[6]

2.4 Contact Force Regulation

Theory

The controller maintains a target contact force $F^*$ along a specified direction $\hat{n}$ through a proxy-motion strategy: the EE position is nudged toward or away from the surface to reduce force error.

Force error: $$e_F = F^* - F_{contact}$$

Compliant approach — binary gate:

  • If $F_{contact} &gt; F_{max}$: stop immediately, return False
  • Otherwise: step EE by $v_{approach} \cdot \Delta t$ along $\hat{n}$

Contact-force maintenance loop (maintain_contact): The position correction magnitude is proportional to $e_F$; direction is inferred from the sign of $e_F$ (too little force → press harder; too much → back off).

Jacobian transpose (conceptual, apply_constant_force): $$\tau = J(\mathbf{q})^T,\mathbf{F}$$

approximated in the current implementation by direct Cartesian EE displacement.

Parameters

Name Default
target_force 5.0 N
force_tolerance 1.0 N
contact_threshold 0.1 N
approach_speed 0.01 m/s

Implementation

fullstack_manip/control/high_level/force_control_controller.py
  ForceControlController(BaseController)
    execute(target_pos, target_orient, force_direction)
    apply_constant_force(force_vector, duration)
    maintain_contact(contact_body, duration) → bool
    compliant_motion(target_pos, max_force, approach_speed=0.01) → bool

2.5 Gravity & Friction Compensation

Theory

Gravity compensation — MuJoCo computes $\boldsymbol{\tau}_{bias}$ containing gravity, Coriolis, and centrifugal torques in data.qfrc_bias. Applying the negative of this term to the actuators makes the robot feel weightless:

$$\boldsymbol{\tau}_{cmd} = k_{comp} \cdot \mathbf{data.qfrc_bias}[:n_q]$$

with $k_{comp} \in [0.0,, 1.5]$ (> 1 slightly over-compensates).

Friction compensation — viscous + Coulomb model:

$$\tau_j^{fric} = -(B_j,\dot{q}_j + \text{sign}(\dot{q}_j) \cdot B_j \cdot 0.5)$$

where $B_j$ is the viscous friction coefficient for joint $j$. Combined command:

$$\boldsymbol{\tau}_{cmd} = \boldsymbol{\tau}_{grav} + \boldsymbol{\tau}_{fric}$$

Friction calibration (calibrate_friction) — collects $(\tau_{applied},, \dot{q})$ pairs at a slow sweep, then estimates:

$$B_j \approx \overline{\left|\frac{\tau_{applied,j}}{|\dot{q}_j|}\right|}$$

Position-holding superimposes a joint-space spring:

$$\boldsymbol{\tau}_{hold} = k_{stiff} \cdot 100 \cdot (q^* - q_{current})$$ $$\boldsymbol{\tau}_{cmd} = \boldsymbol{\tau}_{grav} + \boldsymbol{\tau}_{hold}$$

Parameters

Name Default Range
compensation_gain $k_{comp}$ 1.0 [0.0, 1.5]
enable_friction_compensation False bool
friction_coefficients $B_j$ 0.1 (uniform) per-joint, calibrated

Implementation

fullstack_manip/control/high_level/gravity_compensation_controller.py
  GravityCompensationController(BaseController)
    execute(duration, record_trajectory=False) → [q]
    hold_position(target_joints, duration, stiffness=0.1)
    compensate_and_teach(duration) → [(qpos, ee_pose)]
    calibrate_friction(duration, velocity_amplitude)
    set_compensation_gain(gain)
    _apply_compensation_torques(torques)   # writes data.ctrl, calls mj_step

2.6 Pick-and-Place Sequencer

Theory

A finite state machine with four states driven by Cartesian waypoints and gripper commands. The gripper grasp condition is evaluated by contact force thresholds.

State sequence (pick_object):

# State Action
1 Pre-grasp Move EE to object_pos + [0.04, 0, 0.05] (above)
2 Open Open gripper, verify release
3 Grasp Move EE to object_pos + [0.02, 0, -0.005] (insert)
4 Close Close gripper; if grasp_success == False → retry from step 1

Gripper success is determined by Gripper.check_grasp_success(): contact force between jaw geoms and object geom exceeds grasp_threshold (1.0 N).

Parameters

Offset Value
Pre-grasp height [0.04, 0, 0.05] m
Grasp insert [0.02, 0, -0.005] m
Grasp force threshold 1.0 N
Release force threshold 0.01 N

Implementation

fullstack_manip/control/high_level/pick_place_controller.py
  PickPlaceController(BaseController)
    execute(pick_position, place_position)
    pick_object(object_position) → bool
    place_object(target_position)

Part III — Collision Detection & Avoidance

This part covers all algorithms that detect, measure, and avoid collisions: the shared CollisionCheckerInterface constraint layer used inside the IK QP, three interchangeable geometry backends, and the trajectory-repair gradient-repulsion planner.


3.1 IK Constraint Stacking (LimitManager)

Theory

Mink stacks multiple constraint layers in the QP:

$$\min_{\dot{q}} \tfrac{1}{2}|\dot{q}|^2 \quad \text{s.t.} \quad c_i(\dot{q}) \leq 0 \quad \forall, i \in {\text{config},, \text{collision},, \text{velocity}}$$

Three constraint types are used:

  1. Configuration limit — joint position bounds from URDF: $$q_{lo} \leq q + \dot{q},\Delta t \leq q_{hi}$$

  2. Collision avoidance limit — minimum-distance repulsion between gripper geometry and obstacle geoms. Each pair $(g_{robot}, g_{obs})$ generates an inequality keeping signed distance $d &gt; d_{margin}$: $$\nabla_q d \cdot \dot{q} \geq -\alpha,(d - d_{margin})$$

  3. Velocity limit — per-joint speed cap: $$|\dot{q}j| \leq v{j,max}$$

Default velocity limits (all joints): $v_{max} = \pi$ rad/s.

Implementation

fullstack_manip/core/limit.py
  LimitManager
    set_limits(collision_pairs, max_velocities) → List[mink.Limit]
    get_limits() → List[mink.Limit]

Collision pairs are constructed as:

for body in gripper_bodies:
    for geom_id in mink.get_body_geom_ids(model, body_id):
        for obs_geom_id in obstacle_geom_ids:
            geom_pairs.append(([geom_id], [obs_geom_id]))

3.2 MuJoCo Narrow-Phase Backend

Theory

MuJoCo runs a narrow-phase contact solver during mj_step / mj_forward. It fills data.contact[0:ncon] with contact structs containing:

  • geom1, geom2 — colliding geom IDs
  • dist — signed penetration distance (negative = interpenetrating)
  • frame[:3] — contact normal force resultant (3-vector)

Queries are $O(n_{contacts})$ lookups — no geometry recomputation.

Contact force magnitude:

$$F_{contact} = \sum_{c \in \text{contacts}(A,B)} |\mathbf{c.frame}[:3]|_2$$

Trajectory check — sets data.qpos and calls mj_forward at each waypoint, then scans the contact list for robot–obstacle pairs.

Implementation

fullstack_manip/core/collision_backends/mujoco_collision.py
  MuJoCoCollisionChecker(model, data)
    detect_contact(body_a, body_b) → bool
    check_contact_pair(body_a, body_b) → bool   # dist < 0
    compute_contact_force(body_a, body_b) → float
    minimum_distance(body_a, body_b) → float    # +inf if no contact
    check_trajectory(trajectory[T,nq], obstacle_bodies, robot_bodies=None) → bool
    check_collision(trajectory, obstacles) → bool  # deprecated alias

Internal: _body_geom_ids(name) uses model.body_geomadr[id] + body_geomnum[id].


3.3 GJK/EPA Backend (Coal)

Theory

Coal (fork of hpp-fcl) implements the Gilbert-Johnson-Keerthi (GJK) algorithm for convex objects and the Expanding Polytope Algorithm (EPA) for penetration depth.

GJK iteratively finds the minimum-norm point of the Minkowski difference $A \ominus B$. If the origin lies inside the Minkowski difference, the objects overlap.

Distance query: $$d = \min_{\mathbf{a}\in A,, \mathbf{b}\in B} |\mathbf{a} - \mathbf{b}|$$

BVH mesh representation: BVHModelOBBRSS decomposes a triangle mesh into a hierarchy of oriented bounding boxes (OBB) for $O(\log n)$ broadphase rejection.

FK is provided by yourdfpy: update_cfg(joint_dict)get_transform(link, collision_geometry=True)$T \in SE(3)$. The transform is applied to a Coal object via Transform3s.setRotation(R) + setTranslation(t).

Implementation

fullstack_manip/core/collision_backends/coal_collision.py
  CoalCollisionChecker(urdf_path, joint_names=None)
    detect_contact(body_a, body_b) → bool
    check_contact_pair(body_a, body_b) → bool
    minimum_distance(body_a, body_b) → float
    compute_contact_force(...) → NotImplementedError
    check_trajectory(trajectory, obstacle_bodies, robot_bodies=None) → bool

Install: pip install "fullstack-manip[coal]" (pulls coal>=2.0.0, yourdfpy>=0.0.53).


3.4 BVH Backend (Trimesh)

Theory

Trimesh wraps python-fcl in a high-level CollisionManager. The manager maintains an internal BVH (Bounding Volume Hierarchy) of triangle meshes and supports:

  • in_collision_single(mesh) — broadphase AABB then narrowphase triangle-triangle test
  • min_distance_single(mesh) — minimum unsigned distance between mesh surfaces

Signed distance is reconstructed by negating the unsigned distance when in_collision_single returns True.

FK (same as Coal backend): yourdfpy.update_cfgget_transformmesh.apply_transform(T).

Implementation

fullstack_manip/core/collision_backends/trimesh_collision.py
  TrimeshCollisionChecker(urdf_path, joint_names=None)
    detect_contact(body_a, body_b) → bool
    check_contact_pair(body_a, body_b) → bool
    minimum_distance(body_a, body_b) → float   # negative if penetrating
    compute_contact_force(...) → NotImplementedError
    check_trajectory(trajectory, obstacle_bodies, robot_bodies=None) → bool

Install: pip install "fullstack-manip[trimesh]" (pulls trimesh>=4.0.0, yourdfpy>=0.0.53).


3.5 Gradient Repulsion in C-Space

Theory

The planner repairs a naive joint trajectory $Q = {q^{(i)}}_{i=1}^N$ by pushing colliding waypoints away from the obstacle using a numerical Jacobian-based gradient in configuration space.

Bounding-sphere collision test — fast broadphase per link $L$:

$$\mathbf{c}_L = R_L(\mathbf{q}),\mathbf{c}_{L,local} + \mathbf{t}_L(\mathbf{q})$$ $$d_L = |\mathbf{c}_L - \mathbf{c}_{obs}| + \epsilon_{tiny}$$

Link $L$ collides if $d_L &lt; r_{obs} + r_L$ (sum of bounding-sphere radii).

Numerical Jacobian column for joint $j$ on link $L$:

$$\mathbf{q}^{(j+)} = \mathbf{q},\quad q^{(j+)}_j \mathrel{+}= \varepsilon$$ $$\mathbf{J}^j_L = \frac{\mathbf{c}_L(\mathbf{q}^{(j+)}) - \mathbf{c}_L(\mathbf{q})}{\varepsilon}$$

Repulsion gradient (accumulated over all colliding links):

$$\hat{\mathbf{d}}_L = \frac{\mathbf{c}_L - \mathbf{c}_{obs}}{d_L}, \qquad p_L = (r_{obs} + r_L) - d_L$$ $$\nabla_j \mathrel{+}= (\mathbf{J}^j_L \cdot \hat{\mathbf{d}}_L),p_L$$

Normalised gradient step with joint-limit clamp:

$$\mathbf{q} \leftarrow \mathbf{q} + \alpha,\frac{\nabla}{|\nabla|_2}$$ $$\mathbf{q} \leftarrow \text{clip}(\mathbf{q},, \mathbf{q}_{lo},, \mathbf{q}_{hi})$$

Repeat until collision-free or max_iter steps.

Post-processing:

  1. Gaussian smoothing — removes per-waypoint kinks introduced by local gradient steps. Per-joint convolution with kernel: $$k_i = \frac{\exp(-i^2 / 2\sigma^2)}{\sum_j \exp(-j^2/2\sigma^2)}, \qquad \sigma = 3\ \text{waypoints}$$

  2. Upsample — linear interpolation from planning resolution ($N_{plan}$ waypoints) to playback resolution ($N_{play} = fps \times T$) via np.interp.

  3. Pre-bake — all FK transforms and per-link collision statuses are computed offline. The 60 fps render loop does only $O(n_{links})$ handle-property writes.

Parameters

Parameter Default Description
max_iter 25 Max gradient steps per waypoint
step $\alpha$ 0.05 Gradient step size
eps $\varepsilon$ 0.003 rad Finite-difference perturbation
sigma 3.0 waypoints Gaussian smoothing width
obs_eff_r $|\mathbf{s}/2|_2 + 0.015$ m Obstacle sphere radius + clearance
plan_steps 150 Planning waypoints
fps 60 Playback frame-rate

Complexity

Per waypoint per iteration: $O(n_{links} \times n_{joints})$ FK calls. With plan_steps=150, max_iter=25, n_links=7, n_joints=6, this is at most $150 \times 25 \times 7 \times 6 = 157,500$ FK evaluations (~0.19 ms/call → ≤ 30 s worst-case; typical: 5–8 s, since most waypoints converge in < 5 iterations).

Implementation

examples/collision_vis_viser.py
  _plan_avoiding_trajectory(raw_traj, robot_urdf, jnames,
                             local_centroids, link_radii,
                             obs_center, obs_eff_r,
                             j_lower, j_upper,
                             max_iter=25, step=0.05, eps=0.003)
  _gauss_smooth(traj, sigma=3.0) → smoothed_traj
  _resample(traj, n_out) → resampled_traj
  _hitting_links(local_centroids, link_radii, fk, obs_center, obs_eff_r) → set
  _bake(traj, ...) → List[{q, fk, hits}]

3.6 Backend Comparison

Algorithm Backend Geometry Force data Speed Dependency
MuJoCo narrow-phase MuJoCoCollisionChecker Implicit (sim) Fastest (in-step) mujoco (core)
GJK/EPA CoalCollisionChecker BVH mesh + primitives Fast coal, yourdfpy
BVH trimesh TrimeshCollisionChecker BVH mesh Medium trimesh, yourdfpy
Gradient repulsion _plan_avoiding_trajectory Spheres (broadphase) Planning-time none extra

Select with create_collision_checker(backend, **kwargs) from fullstack_manip.core.collision.


Part IV — Trajectory Generation & Optimization

This part covers algorithms in fullstack_manip/planning/trajectory_generators/ and planning/task/.
Sections marked [🔲 not yet implemented] are still design specifications; see
sota-planning.md for library recommendations.


4.1 Minimum-Jerk Polynomial

Theory

A degree-5 polynomial minimises the integral of squared jerk:

$$\min \int_0^T \dddot{x}^2,dt \quad \Rightarrow \quad x(t) = x_0 + (x_f - x_0)\left[10\tau^3 - 15\tau^4 + 6\tau^5\right], \quad \tau = \frac{t}{T}$$

Velocity and acceleration are guaranteed zero at both endpoints. Generalises to joint-space by applying the formula per-DOF.

Flash & Hogan (1985) showed this polynomial matches biological reaching trajectories.

Parameters

Name Value Description
T caller-specified Duration in seconds
n_joints 6 Applied per-DOF

Implementation

fullstack_manip/planning/trajectory_generators/min_jerk.py
  min_jerk_trajectory(q_start, q_end, duration, dt=0.01) → (traj[T,n], t[T])
  min_jerk_path(waypoints, durations?, dt=0.01, default_speed=0.5) → (traj, t)
  MinJerkGenerator(dt=0.01, default_speed=0.5)
    .generate(q_start, q_end, duration?) → (traj, t)
    .generate_path(waypoints, durations?) → (traj, t)
from fullstack_manip.planning import min_jerk_trajectory, min_jerk_path, MinJerkGenerator

# Single segment
traj, t = min_jerk_trajectory(q0, q1, duration=1.0)

# Piecewise multi-waypoint (duration auto-estimated from joint-space distance)
traj, t = min_jerk_path([q0, q_via, q1], default_speed=0.5)

# Stateful wrapper
gen = MinJerkGenerator(dt=0.01, default_speed=0.5)
traj, t = gen.generate(q0, q1)           # duration = ‖q1-q0‖₂ / default_speed
traj, t = gen.generate_path([q0, q1])    # piecewise

4.2 Time-Optimal Parametrisation (TOPP-RA)

Theory

Given a geometric path $q(s), s \in [0, 1]$ (e.g. from a cubic spline), find the
fastest time parametrisation $s(t)$ respecting per-joint velocity and acceleration limits:

$$\min_{\dot{s}\geq 0} \int_0^1 \frac{ds}{\dot{s}} \quad \text{s.t.}\quad |\dot{q}(s)\dot{s}| \leq \dot{q}^{\max},\quad |\ddot{q}(s)\dot{s}^2 + q'(s)\ddot{s}| \leq \ddot{q}^{\max}$$

TOPP-RA (Pham & Pham, 2018) solves this via a 1-D reachability analysis scan — $O(N)$ in path waypoints.

Library

pip install toppra

Implementation

fullstack_manip/planning/trajectory_generators/topp_ra.py
  reparametrize(trajectory, t_in, velocity_limits?, acceleration_limits?, gridpoints?)
      → (new_traj, new_t)
  ToppRaGenerator(velocity_limits?, acceleration_limits?)
    .reparametrize(trajectory, t_in) → (new_traj, new_t)
    .available → bool   # False when toppra not installed → returns (traj, t) unchanged
from fullstack_manip.planning import ToppRaGenerator, reparametrize

new_traj, new_t = reparametrize(traj, t, velocity_limits=np.full(6, np.pi))

gen = ToppRaGenerator(velocity_limits=np.full(6, np.pi),
                      acceleration_limits=np.full(6, 2*np.pi))
new_traj, new_t = gen.reparametrize(traj, t)

Graceful degradation: emits UserWarning at import when toppra is missing; call-site code unchanged.
Integrates with: any trajectory as a post-processing step.


4.3 STOMP Trajectory Optimization

Theory

STOMP (Stochastic Trajectory Optimization for Motion Planning, Kalakrishnan et al. 2011) generates noisy trajectory samples, evaluates cost (smoothness + obstacle penalty), and updates the trajectory via importance-weighted average:

$$\theta_{new} = \theta + \frac{\sum_k \exp(-\frac{1}{\lambda}S(\Theta_k)) \cdot \delta\Theta_k}{\sum_k \exp(-\frac{1}{\lambda}S(\Theta_k))}$$

where $\Theta_k = \theta + \epsilon_k$ are perturbed samples and $S$ is the total cost.

No gradient computation required — purely sample-based.

Parameters

Name Default Description
n_samples 20 Number of noisy trajectories per iteration
n_iter 50 Optimization iterations
noise_stddev 0.05 rad Per-DOF Gaussian noise
lambda 10.0 Temperature controlling exploitation
obstacle_weight 1.0 Weight of obstacle penetration cost

Implementation

[🔲 not yet implemented] — target file: fullstack_manip/planning/trajectory_optimizers/stomp.py

# Planned API
def stomp_optimize(trajectory, is_state_valid_fn, n_samples=20, n_iter=50) -> np.ndarray:
    """Returns optimized trajectory [T, n_joints]."""
    ...

4.4 iLQR / Differential Dynamic Programming

Theory

iLQR (iterative Linear Quadratic Regulator) optimises a trajectory under nonlinear dynamics
by alternating backward-pass (compute linear-quadratic gains) and forward-pass (improve trajectory):

$$\min_{x_{1:T}, u_{1:T}} \sum_{t=1}^{T-1} l(x_t, u_t) + l_f(x_T)$$

The backward pass solves a Riccati recursion; the forward pass does a line search along the optimal perturbation direction. DDP adds second-order dynamics terms for faster convergence.

Library

pip install crocoddyl  # or implement from scratch with JAX/NumPy

Implementation

[🔲 not yet implemented] — suitable for torque-controlled arms in MuJoCo simulation;
not directly applicable to position-mode soarm100 servos.

Target integration: OperationalSpaceController torque-mode (already scaffolded) →
iLQR provides the reference trajectory; OSC tracks it.


4.5 Ground-Atom STRIPS Task Planner

Theory

A ground-atom STRIPS planner represents task state as a set of Boolean atoms (strings). Each TaskAction specifies:

  • preconditions — atoms that must all hold before the action can fire
  • effects_add — atoms that become true after the action
  • effects_del — atoms that become false after the action

Search uses breadth-first search (BFS) over the reachable state space, guaranteeing the shortest plan (fewest actions) when one exists within max_depth.

$$\text{plan}^* = \arg\min_{\pi} |\pi| \quad \text{s.t.} \quad \text{goal} \subseteq \text{result}(\text{init},, \pi)$$

Complexity: $O(b^d)$ where $b$ is the branching factor (number of applicable actions at each state) and $d$ is the plan depth. Practical for manipulation sequences (typically $d \leq 10$, $b \leq 6$).

Parameters

Name Default Description
actions MANIPULATION_OPERATORS List of TaskAction objects
max_depth 50 BFS depth cap — prevents unbounded search

Predefined Operators (MANIPULATION_OPERATORS)

Action Preconditions Effects
open_gripper gripper_closed +gripper_open, -gripper_closed
close_gripper gripper_open, arm_at_grasp_pose +gripper_closed, +object_grasped, -gripper_open
move_to_grasp_pose gripper_open +arm_at_grasp_pose
move_to_place_pose object_grasped +arm_at_place_pose
move_to_home +arm_at_home
release_object gripper_closed, arm_at_place_pose +gripper_open, +object_placed, -object_grasped

Implementation

fullstack_manip/planning/task/task_planner.py
  TaskAction(name, preconditions, effects_add, effects_del, metadata)
    .is_applicable(state) → bool
    .apply(state) → new_state
  TaskPlanner(actions?, max_depth=50)
    .plan(initial_state, goal_state) → TaskPlan
    .add_action(action)   # replaces by name
  MANIPULATION_OPERATORS  — predefined tabletop pick-place ground actions

fullstack_manip/core/interfaces.py
  TaskPlan(success, actions, info)   — result dataclass
  TaskPlannerInterface               — runtime-checkable Protocol
from fullstack_manip.planning import TaskPlanner, MANIPULATION_OPERATORS

planner = TaskPlanner(MANIPULATION_OPERATORS, max_depth=10)
plan = planner.plan(
    initial_state={'gripper_open'},
    goal_state={'object_placed'},
)
print(plan.success)   # True
print(plan.actions)   # ['move_to_grasp_pose', 'close_gripper',
                      #  'move_to_place_pose', 'release_object']

Part V — Analytics & Telemetry

fullstack_manip/analysis/ provides opt-in, zero-overhead telemetry instrumentation for controllers and planners. Every instrumented component accepts a telemetry_logger=None argument; omitting it disables all recording with no runtime cost.


5.1 TelemetryBuffer

Theory

A thread-safe, in-memory ring buffer that stores (timestamp, key, value) triples. The buffer is backed by collections.deque(maxlen=N) so older entries are automatically discarded when the buffer is full, making it suitable for long-running experiments.

Data Model

TelemetryEntry = (timestamp: float, key: str, value: Any)
TelemetryBuffer.data: deque[TelemetryEntry]

API

from fullstack_manip.analysis import TelemetryBuffer

buf = TelemetryBuffer(maxlen=5000)   # ring buffer, oldest entries dropped when full

# Log a scalar / array / dict
buf.log("joint_pos", q)              # (6,) array
buf.log("ee_force", wrench)
buf.log("step", {"q": q, "dq": dq})

# Retrieve as NumPy array
q_history  = buf.to_numpy("joint_pos")   # shape (N, 6)
timestamps = buf.timestamps("joint_pos") # shape (N,)

# Persist
buf.save("run_001.npz")
buf2 = TelemetryBuffer.load("run_001.npz")

# Introspection
len(buf)                             # total entries across all keys
buf.keys()                           # set of logged key names
buf.clear()

Implementation Notes

  • log() acquires a threading.Lock — safe to call from control thread
  • to_numpy(key) stacks all values matching key; raises KeyError if key absent
  • save() persists to .npz with three arrays per key: timestamps_<key>, values_<key>, and a _keys metadata array
  • load() reconstructs the deque in chronological order

5.2 Telemetry Logger Protocol

Protocol Definition

# fullstack_manip/core/interfaces.py

class TelemetryLoggerInterface(Protocol):
    def log(self, key: str, value: Any, timestamp: float | None = None) -> None: ...

Any class implementing log(key, value, timestamp) satisfies the protocol. No base class required; duck typing applies.

NullLogger

from fullstack_manip.analysis import NullLogger

logger = NullLogger()
logger.log("anything", 42)   # silent no-op

NullLogger is the drop-in when telemetry is disabled. It carries zero overhead because Python's method dispatch is not guarded by an if logger is not None check in hot paths — the _emit() helper on BaseController handles this transparently (see §5.5).


5.3 RerunLogger

Theory

RerunLogger wraps the Rerun SDK to stream telemetry in real time to the Rerun desktop viewer. Useful for visualising 3-D trajectories, joint time-series, and contact forces during live experiments.

Entity Path Convention

Key logged Rerun entity path Type logged
joint_pos robot/joints/position rr.Scalar per joint
joint_vel robot/joints/velocity rr.Scalar per joint
ee_pos robot/ee/position rr.Points3D
ee_force robot/ee/force rr.Scalar
rrt_node planning/rrt/nodes rr.Points3D
trajectory planning/trajectory rr.LineStrips3D
controller_error control/error rr.Scalar

Custom keys are forwarded to <key> entity path as rr.AnyValues.

API

from fullstack_manip.analysis import RerunLogger

logger = RerunLogger(
    application_id="arm_debug",         # Rerun viewer window title
    recording_id=None,                  # auto-generated UUID if None
    spawn=True,                         # spawn Rerun viewer process inline
    blueprint=None,                     # optional rr.Blueprint for layout
)

# Attach to a component
planner = RRTPlanner(..., telemetry_logger=logger)

Graceful Degradation

If rerun-sdk is not installed the module-level _HAS_RERUN guard disables the logger silently:

try:
    import rerun as rr
    _HAS_RERUN = True
except ImportError:
    _HAS_RERUN = False

RerunLogger.log() becomes a no-op when _HAS_RERUN is False.


5.4 Offline Plot Functions

All plot functions live in fullstack_manip/analysis/ and return a matplotlib.figure.Figure. They do not call plt.show(), so figures can be saved, embedded in notebooks, or forwarded to the Streamlit dashboard without side-effects.

plot_trajectory.py

Function Input keys Description
plot_joint_positions(buf_or_array) joint_pos (N,6) Joint angles vs time, one axis per joint
plot_joint_velocities(buf_or_array) joint_vel (N,6) Joint velocities vs time
plot_ee_path(buf_or_array) ee_pos (N,3) 3-D end-effector path (Axes3D)
plot_trajectory_overview(buf) joint_pos, ee_pos Combined 2×2 summary figure

plot_controller.py

Function Input keys Description
plot_command_vs_actual(buf) cmd, actual Commanded vs measured positions
plot_force_signals(buf) ee_force Wrench components vs time
plot_controller_error(buf) controller_error Tracking error norm vs time
plot_controller_overview(buf) all above Combined summary figure

plot_planning.py

Function Input keys Description
plot_rrt_tree(buf) rrt_node Scatter of explored C-space nodes
plot_cost_convergence(buf) cost Cost / path length vs iteration
plot_planning_overview(buf) all above Combined planning summary figure

Usage

from fullstack_manip.analysis import TelemetryBuffer
from fullstack_manip.analysis.plot_trajectory import plot_joint_positions, plot_ee_path
from fullstack_manip.analysis.plot_controller import plot_controller_overview
from fullstack_manip.analysis.plot_planning import plot_planning_overview

buf = TelemetryBuffer.load("run_001.npz")

fig1 = plot_joint_positions(buf)
fig1.savefig("joints.png", dpi=150)

fig2 = plot_ee_path(buf)
fig2.savefig("ee_path.png", dpi=150)

5.5 Controller Telemetry Emission Pattern

_emit() Helper

BaseController adds the _emit() helper that routes through the injected logger:

# fullstack_manip/control/high_level/base_controller.py (simplified)

class BaseController:
    def __init__(self, ..., telemetry_logger=None):
        self._tel = telemetry_logger  # None → no-op

    def _emit(self, key: str, value) -> None:
        if self._tel is not None:
            self._tel.log(key, value)

Every controller subclass calls self._emit() inside execute():

# ImpedanceController
def execute(self, target_pos, target_rot=None, **kwargs):
    ...
    self._emit("joint_pos", q)
    self._emit("ee_pos", ee_pos)
    self._emit("controller_error", error_norm)
    ...

Instrumented Components

Component Keys emitted
ImpedanceController joint_pos, ee_pos, controller_error
TrajectoryFollowingController joint_pos, joint_vel, cmd
MotionPlanner.plan_trajectory trajectory, joint_pos
RRTPlanner.plan rrt_node, cost

Injecting a Logger into a Controller

from fullstack_manip.analysis import TelemetryBuffer
from fullstack_manip.control.high_level.impedance_controller import ImpedanceController

buf = TelemetryBuffer(maxlen=10_000)
ctrl = ImpedanceController(robot=robot, dt=0.01, telemetry_logger=buf)

for _ in range(1000):
    ctrl.execute(target_pos=goal)

# Post-process
from fullstack_manip.analysis.plot_controller import plot_controller_overview
fig = plot_controller_overview(buf)
fig.savefig("ctrl_analysis.png")