Certifiably Optimal LEO Doppler Positioning

• The paper introduces a certifiably optimal estimation method that transforms the nonlinear LEO Doppler positioning problem into a convex semidefinite programming framework.
• It employs graduated weight approximation and QCQP reformulation to eliminate local minima and ensure global optimality even with low-noise measurements.
• Simulations and real-world tests demonstrate that the SDP-based approach outperforms traditional methods, achieving 3D errors as low as 140 m and providing reliable initialization-free positioning.

A certifiably optimal Low Earth Orbit (LEO) Doppler positioning method is an estimation algorithm that produces globally optimal user position and clock drift estimates from Doppler shift measurements collected from LEO satellites, with mathematical guarantees for global optimality under controllable noise and modeling conditions. This approach systematically eliminates the risk of convergence to local minima inherent to conventional iterative solvers and is designed for robust, initialization-free positioning, even with signals of opportunity (“SOPs”) in unknown environments.

1. Problem Formulation and Model Structure

The fundamental measurement model expresses the normalized Doppler observation $D_i$ (for satellite $i$) as: $$D_i = \frac{(\mathbf{p}_r - \mathbf{p}^s_i)^\top (\mathbf{v}_r - \mathbf{v}^s_i)}{\rho_i} + c \cdot \dot{d}_t + \epsilon$$ where:

• $\mathbf{p}_r$ and $\mathbf{v}_r$ are the receiver’s unknown position and velocity,
• $\mathbf{p}^s_i$, $\mathbf{v}^s_i$ are the known position and velocity of satellite $i$,
• $\rho_i = \|\mathbf{p}_r - \mathbf{p}^s_i\|$ is the geometric range,
• $c$ is the speed of light,
• $\dot{d}_t$ is the receiver clock drift, and
• $\epsilon$ is the additive measurement noise (often Gaussian).

The aim is to estimate $\mathbf{p}_r$ and $\dot{d}_t$ from a set of $N$ such measurements, minimizing the nonlinear weighted least-squares (NWLS) objective: $$\min_{\mathbf{p}_r,\, \dot{d}_t} \sum_{i=1}^N \left\| D_i - \frac{(\mathbf{p}_r - \mathbf{p}^s_i)^\top (\mathbf{v}_r - \mathbf{v}^s_i)}{\rho_i} - c \dot{d}_t \right\|^2_Q$$ where $Q$ is a weighting matrix, generally diagonal with entries $1/\rho_i$.

This problem is nonconvex, due to the fractional term and the geometric range dependence on $\mathbf{p}_r$, and has multiple local minima (Song et al., 21 Sep 2025).

2. Convexification via Graduated Weight Approximation (GWA) and Semidefinite Relaxation

2.1 Fractional to Polynomial Reformulation

The original NWLS objective is multiplied by the denominators $\rho_i$ to yield a polynomial optimization problem (POP), facilitating treatment with polynomial optimization techniques. The weights are updated iteratively according to the latest estimated ranges in a process called graduated weight approximation (GWA).

Let $Q$ be frozen for each GWA iteration, treating the problem as: $$\min_{\mathbf{p}_r,\, \dot{d}_t} \sum_{i=1}^N \left( \rho_i D_i - (\mathbf{p}_r - \mathbf{p}^s_i)^\top (\mathbf{v}_r - \mathbf{v}^s_i) - c \rho_i \dot{d}_t \right)^2$$ where monomials in $\mathbf{p}_r$ and $\rho_i$ can be expanded as polynomials.

2.2 Lifting and Quadratically Constrained Quadratic Program (QCQP) Formulation

Introduce auxiliary variables $\rho_i$ and $z_i = c \,\rho_i \,\dot{d}_t$, with explicit algebraic constraints linking these variables. The decision variable becomes: $$\mathbf{y} = [\mathbf{p}_r^\top,\, c\dot{d}_t,\, \rho_1,\, \ldots,\, \rho_N,\, z_1,\, \ldots,\, z_N]^\top$$ With these, the objective and constraints can be written as quadratic forms or quadratic constraints, yielding a QCQP that is still nonconvex (due to the quadratic equality $Y = \mathbf{y}\mathbf{y}^\top$ implicit in the lifting).

2.3 Semidefinite Programming (SDP) Relaxation

The key step is relaxing the nonconvex $Y = \mathbf{y}\mathbf{y}^\top$ to $Y \succeq \mathbf{y}\mathbf{y}^\top$ and formulating a linear matrix inequality: $$S=\begin{bmatrix} Y & \mathbf{y} \ \mathbf{y}^\top & 1 \end{bmatrix} \succeq 0$$ while dropping the rank-1 constraint on $Y$—the standard Shor SDP relaxation. The resulting problem is convex and solvable to global optima with polynomial-time SDP solvers.

3. Theoretical Optimality Guarantees

3.1 Noiseless Case

In the noise-free case, necessary optimality conditions (Karush-Kuhn-Tucker, KKT) for the original NWLS are:

• All constraint functions vanish at the optimum,
• The dual matrix $H(\lambda)$ is positive semidefinite,
• The complementarity condition $H(\lambda)[1; \mathbf{y}^*]=0$ is met,
• The rank of $H(\lambda)$ is $n-1$ (where $n$ is the number of variables in $\mathbf{y}$).

These conditions guarantee that the SDP solution is tight: the optimizer matrix $Y$ has rank-1 and can be factored as $Y = \mathbf{y}^*\mathbf{y}^{*\top}$. Thus, the global minimizer is uniquely recovered from the SDP solution.

3.2 Noisy Case

Global optimality is preserved under sufficient control of the noise level. Specifically, if the measurement noise is lower than a threshold specified in terms of the singular values of the Jacobian (of the constraints) and the spectral gap of the quadratic form matrix: $$\frac{1}{\sigma_N}\|\mathcal{G}\|\, \|\nabla f_\theta(\bar{y})\| + \|F-F_\theta\| < \nu_{N+4}(F_\theta)$$ then the SDP relaxation remains tight, guaranteeing that a rank-1 solution is recovered [(Song et al., 21 Sep 2025), cf. Proposition 1].

4. Simulation and Experimental Benchmarking

4.1 Monte Carlo Simulations

The SDP relaxation-based method is benchmarked against Gauss-Newton (GN) and Dog-Leg (DL) local solvers across a range of initialization distances $\Delta_0$.

• For small $\Delta_0$ ($<$100 km), all methods converge to the global optimum.
• As $\Delta_0$ increases ($\gtrsim$1000 km), GN and DL regularly converge to spurious local optima or fail, with position errors in the range of 1,000 to 2,000 km.
• The SDP always returns the global optimum, independent of initialization, with 3D error $\approx$0.71 km.

Further refinement by using the SDP solution as the initializer for GN/DL yields $\approx$0.13 km 3D error, demonstrating that tight convexification not only certifies optimality but also improves subsequent local refinement.

4.2 Real-World Data

On a real test path using 8 Iridium-NEXT satellites (35 seconds observation), local methods (GN, DL) with $\Delta_0>$1000 km produce solutions thousands of kilometers off the ground truth, while the SDP-based solver consistently yields a certifiably global solution with 3D error $\approx$140 m. Using the SDP output as initialization for local search lowers the error to $\approx$130 m.

5. Implementation Considerations and Practical Impact

• Initialization-Free: Given its convexity, the SDP method removes the need for an accurate initial guess—a critical property for navigation in environments where the user state is a priori unknown.
• Certifiability: A posteriori, global optimality is guaranteed by inspection of the SDP duality gap and matrix rank, delivering a “certificate” of optimality.
• Integration with Local Search: Under mild noise, the SDP solution is both optimal and computationally tractable; under higher noise or when strict tightness cannot be confirmed, the SDP solution provides an excellent initialization for local optimization, which typically brings further refinement.
• Noise Bound Verification: The noise bounds are explicit; if measurement noise exceeds the bound required for tightness, global optimality is not guaranteed and diagnostics will indicate this, prompting robustification.

6. Key Mathematical Formulas

Concept	Formula/Construction	Description/Notes
Doppler model (single satellite)	$$D=\frac{(\mathbf{p}_r-\mathbf{p}^s)^\top(\mathbf{v}_r-\mathbf{v}^s)}{\rho}+c\dot{d}_t+\epsilon$$	Fundamental measurement model
Range definition	$\rho=\|\mathbf{p}_r-\mathbf{p}^s\|$	Geometric range between receiver and satellite
QCQP constraint (lifting)	$Y=\mathbf{y}\mathbf{y}^\top$ (relaxed to $Y\succeq \mathbf{y}\mathbf{y}^\top$)	Key step for SDP relaxation
SDP feasibility matrix	$$S=\begin{bmatrix} Y&\mathbf{y}\ \mathbf{y}^\top&1\end{bmatrix}\succeq0$$	Forces convexity and enables global solution extraction
Noise bound (for tightness)	$$\frac{1}{\sigma_N}\|\mathcal{G}\| \|\nabla f_\theta(\bar{y})\| + \|F-F_\theta\| < \nu_{N+4}(F_\theta)$$	Certifies that noise is small enough for global recovery

7. Significance and Applications

The certifiably optimal LEO Doppler positioning method outlined here is especially significant for:

• GNSS outage backup: It leverages SOPs from LEO satellites for navigation without requiring strong prior initialization.
• Autonomous navigation in unknown environments: No reliance on a-priori state or environment constraints.
• Network-centric applications: Can be used to cold-start other navigation algorithms, serving as a robust fallback for both civilian and defense platforms.
• Scalability: The SDP approach is naturally parallelizable and can be efficiently solved using state-of-the-art convex solvers.

The approach is validated with both simulation and real-world datasets, confirming that it finds the global solution where conventional solvers fail, and the global solution serves as an effective initializer for rapid, high-precision local refinement (Song et al., 21 Sep 2025).

---

In summary, certifiably optimal LEO Doppler positioning leverages convex optimization—formulated via SDP relaxation after GWA and lifting—to guarantee recovery of the global optimum for user state estimation from Doppler measurements, providing formal mathematical guarantees, practical robustness, and benchmarking that decisively outperforms local search methods in initialization-stressed or ambiguous navigation scenarios.