Q-law GVE & Gradient Validation (#723)

Motivation

Our Q-law implementation (qlaw.py) uses hand-derived Gauss Variational Equation (GVE) coefficients and Q-gradient partials to compute optimal thrust direction. Analytical derivatives can harbour sign errors, missing terms, or incorrect chain-rule application. This document records the three-layer validation performed against independent references.

Q-function Formulation Note

Our Q-function (proximity_quotient) uses a simplified formulation:

Q = Σ W_j × (Δoe_j / S_j)²

with an additive inverse-square periapsis barrier, where S_j = max(|Δoe_j|, floor_j).

This differs from the Petropoulos/pyqlaw formulation which uses multiplicative exponential barriers and oe_dot_max scaling. Consequently, direct Q-gradient comparison with pyqlaw is not meaningful. However, the GVE coefficients are standard orbital mechanics (Battin / Vallado) and must match any correct implementation. Additionally, our analytical gradients can be validated against numerical finite-difference gradients of our own Q-function.

Validation Layers

1. GVE Coefficient Verification

Compare gve_coefficients() against an independently-coded reference implementation of the standard Keplerian GVE (Battin eq. 9.24, Vallado ch. 9):

Element	f_r	f_θ	f_h
da/dt	2a²e·sin(ν)/h	2a²p/(rh)	0
de/dt	p·sin(ν)/h	((p+r)cos(ν)+re)/h	0
di/dt	0	0	r·cos(θ)/h
dΩ/dt	0	0	r·sin(θ)/(h·sin i)
dω/dt	-p·cos(ν)/(he)	(p+r)sin(ν)/(he)	-r·sin(θ)cos(i)/(h·sin i)

The reference implementation recomputes all intermediates from scratch (p, h, r, θ) to avoid sharing any code paths with the production implementation.

Test states (5 representative orbits × 6 true anomalies each = 30 evaluations):

1. LEO circular: a = 6,800 km, e = 0.001, i = 28.5°, Ω = 45°, ω = 0°
2. GTO eccentric: a = 24,400 km, e = 0.73, i = 7°, Ω = 120°, ω = 180°
3. High-inclination: a = 7,200 km, e = 0.01, i = 85°, Ω = 270°, ω = 90°
4. Heliocentric: a = 1.5 AU, e = 0.2, i = 5°, Ω = 100°, ω = 250°
5. Near-equatorial: a = 42,164 km, e = 0.05, i = 0.5°, Ω = 0°, ω = 45°

True anomalies tested: ν ∈ {0°, 45°, 90°, 135°, 180°, 270°}.

Tolerance: < 1e-10 relative error (limited only by float64 representation).

2. Q-Gradient Numerical Validation

Central finite-difference validation of _q_gradient_partials():

∂Q/∂oe_j ≈ (Q(oe_j + ε) − Q(oe_j − ε)) / (2ε)

Each of the 5 test states is perturbed in each of the 5 orbital elements. Step sizes:

• a: ε = 1.0 m (1e-6 × a scale)
• e: ε = 1e-8
• i, Ω, ω: ε = 1e-8 rad

Tolerance: < 1e-6 relative error (limited by ε choice and float cancellation).

3. Periapsis Barrier Gradient Validation

When periapsis radius rp = a(1−e) is between min_rp and rp_safe, the barrier P(rp) = k / (rp − min_rp)² adds gradient contributions via chain rule:

dQ_da += dP/drp × (1 − e)
dQ_de += dP/drp × (−a)

Validated by finite-differencing the full proximity_quotient() + barrier function with barrier-activating states (rp slightly above min_rp).

Tolerance: < 1e-4 relative error (barrier is stiff near min_rp; larger ε needed).

Results

All three validation layers pass at the specified tolerances. No sign errors, missing terms, or chain-rule mistakes found in the production implementation.

Layer	States	Evaluations	Max Relative Error	Status
GVE coefficients	5 × 6ν	450 components	2.60e-16	PASS
Q-gradient partials	5	25 elements	6.55e-09	PASS
Periapsis barrier gradient	1	2 elements (a, e)	2.97e-07	PASS

Note on Q-function scale dependence

The a floor in proximity_quotient is 0.01 × a (depends on a itself). When |error| > floor, the Q contribution saturates at W (constant), making the true derivative 0 — but the analytical gradient gives 2W/Δa because it treats S as fixed. This is correct for steering (the gradient always points toward smaller errors) but means numerical validation must freeze scales at the unperturbed state. Other elements (e, i, Ω, ω) have constant floors and match to full float64 precision.

Validation Artefacts

• Tool: services/tick-engine/tools/qlaw_gradient_check.py — standalone, re-runnable
• Regression tests: TestGradientValidation class in tests/test_qlaw.py

References

• Battin, R.H., An Introduction to the Mathematics and Methods of Astrodynamics, eq. 9.24
• Vallado, D.A., Fundamentals of Astrodynamics and Applications, ch. 9
• Petropoulos, A.E., “Low-thrust orbit transfers using candidate Lyapunov functions with a mechanism for coasting”, AIAA 2004-5089