Pseudorange and Doppler-Based Positioning: Enabling Convergence of Least-Squares Estimation from MEO to LEO

2.1 The LS Problem

In its most general form, solving an LS problem involves addressing the following minimization problem (Dennis & Schnabel, 1996; Quarteroni et al., 2014):

$\underset{x\in {ℝ}^n}{\min }\Phi (x)with\Phi \left(x\right)=\frac{1}{2}{\Vert R(x)\Vert }_2^2=\frac{1}{2}\sum _{i=1}^m{r}_i^2(x)=\frac{1}{2}\sum _{i=1}^m{[y_i-f_i(x\left)\right]}^2$1

where x is a decision variable vector in ℝⁿ, containing the variables to be optimized to minimize the objective function Φ(x). The notation ${\Vert R(x)\Vert }_2^2$ denotes the squared ℓ²-norm of the function $R\left(x\right)={\left(r_1\right(x),\cdots ,r_m(x\left)\right)}^{\top }$ that contains the residuals r_i = y_i – f_i(x). These residuals represent the difference between the observed value y_i and the value predicted by the mathematical model f_i(x) serving as a metric for the discrepancy between the estimated parameters and observed values. The objective of the LS method is to minimize the sum of squares of the residuals, thereby identifying the parameters that most accurately fit the observations. If the function f_i(x) exhibits a linear relationship with the variable x, the problem in Equation (1) is referred to as a linear LS problem; conversely, if there exists a nonlinear relation, the problem will be classified as a nonlinear LS problem. It is important to note that Equation (1) relies on the assumption that the observation noise is independent and identically distributed, an assumption that has been adopted, without a loss of generality, in this work.

In radionavigation and radiolocalization problems as well as in the processing of GNSS observables, the LS is commonly employed to estimate the receiver or target state. In this context, the term x in Equation (1) represents the state to be estimated, which can be written in vector form as follows:

$x=\left[\begin{array}{c}{r}_u\\ v_u\\ {\tau }_u\end{array}\right]$ 2

where r_u denotes a 3 x 1 vector indicating the user’s position relative to a specific reference frame, v_u represents a 3 x 1 vector denoting the user’s velocity, and ${\tau }_u={\left[\delta t_u\delta {\dot{t}}_u\right]}^{{}^{\top }}$ is a 2 x 1 vector containing an unknown local clock offset and drift relative to a specified time scale. Specifically, the clock bias δt_u adjusts for the time offset between the user and the GNSS time scale, whereas the clock drift $\delta {\dot{t}}_u$ determines whether the local oscillator operates at a slower or faster frequency compared with the GNSS system time.

The residuals depend on both the observations and the mathematical model employed in the estimation problem. In this paper, we focus on pseudorange and Doppler shift observations, using the generic observable notation y_i for both cases to ensure uniformity, although the two types of measurements are treated separately throughout the analysis. In the case of pseudorange observations, y_i = ρ_i denotes pseudorange measurements from the i-th satellite, and f_i represents the pseudorange model, expressed as follows:

$f_{{}_i}(x)=\Vert r_i-r_u\Vert +c\delta t_u$3

where r_i is a 3 x 1 vector indicating the position of the i-th satellite and c is the speed of light.

In the case of Doppler shift observations, $y_i={\dot{\rho }}_i$ represents the Doppler shift measurements associated with the signal arriving from the i-th satellite, and f_i denotes the Doppler shift model, mathematically described as follows:

$f_{{}_i}(x)=-\frac{1}{\lambda }\left[\right(v_i-v_u)\cdot \frac{r_i-r_u}{\Vert r_i-r_u\Vert }+c\delta {\dot{t}}_u]$ 4

where v_i is a 3 x 1 vector indicating the velocity of the i-th satellite and λ is the wavelength of the transmitted signal. From Equations (3) and (4), it is evident that a nonlinear relationship exists between the observation and the receiver state. Therefore, our focus shifts to mathematical methods tailored to address such nonlinearity. Whereas linear equations allow direct solutions, the nonlinearity of Φ(x) makes closed-form solutions unfeasible. Therefore, we will investigate iterative techniques, such as Newton’s method, to approximate the solutions.

2.2 Newton Methods to Solve LS Problems

Let us suppose that we want to compute the zeros of a generic function Λ(x). Newton’s method is a numerical technique that leverages the information provided by the function’s derivative to iteratively approximate the zeros. The process begins with an initial guess for the root, which is progressively refined at each iteration through the application of a correction term, δx, to improve the approximation. This correction term is calculated as the ratio of the function Λ(x) to its gradient ∇Λ(x) and is applied to the current estimate x^(k) guiding the algorithm toward a more accurate solution, as shown in the following equation:

$x^{(k+1)}=x^{\left(k\right)}+\delta x^{\left(k\right)}with\delta x^{\left(k\right)}=\frac{\Lambda \left(x^{\left(k\right)}\right)}{\nabla \Lambda \left(x^{\left(k\right)}\right)}$ 5

where the superscript k represents the iteration index. To address the LS problem outlined in Equation (1), it is necessary to examine the stationary points of the objective function Φ(x). This process reduces the problem to solving the equation ∇Φ(x) = 0. It is precisely in looking for the zeros of the gradient of the objective function that Newton’s method proves to be a viable option. When the model in Equation (5) is fitted to our scenario, it becomes evident that Λ(x) = ∇Φ(x) and that ∇Λ(x) = ∇²Φ(x), corresponding to the Hessian of the objective function Φ(x). In light of these considerations, it is possible to reformulate Equation (5) as follows:

• Download figure
• Open in new tab
• Download PowerPoint

where ∇Φ(x^(k)) and ∇²Φ(x^(k)) can be further developed as follows (Gander et al., 2014; Quarteroni et al., 2014):

$\nabla \Phi \left(x^{\left(k\right)}\right)=-J_R{\left(x^{\left(k\right)}\right)}^{\top }R\left(x^{\left(k\right)}\right)=-J_f{\left(x^{\left(k\right)}\right)}^{\top }R\left(x^{\left(k\right)}\right)$7a

${\nabla }^2\Phi \left(x^{\left(k\right)}\right)=J_f{\left(x^{\left(k\right)}\right)}^{\top }{J}_f\left(x^{\left(k\right)}\right)+\sum _{i=1}^m\frac{\partial r_i\left(x^{\left(k\right)}\right)}{\partial x_l{x}_j}{r}_l\left(x^{\left(k\right)}\right),l,j=1,\cdots ,n$7b

Here, the superscript ⊤ denotes the transpose operator. 𝕁_f represents the Jacobian matrix, where each row contains the gradient vector of the corresponding function f_i with respect to the variables of the vector x. Regarding the Hessian matrix in Equation (7b), it is noteworthy that its expression includes two terms. The first term, which is usually the most significant, depends solely on the Jacobian matrix, whereas the second term consists of second-order derivatives multiplied by the residuals. When the model closely approximates the solution, resulting in small residuals, the influence of this second term decreases compared with the first. In addition, computing m × n² derivatives for the second term can be computationally demanding, especially as the dimension of the state increases. Consequently, in practical applications, it is common practice to approximate or even neglect the second term. In particular, the GN method is a well-known method that neglects the second term of Equation (7b), approximating the entire matrix with only the term related to the Jacobian matrix $\left({\nabla }^2\Phi \right(x)\approx J_f{\left(x\right)}^{\top }{J}_f(x\left)\right)$. For this variant, the formulation in Equation (5) can be rewritten as follows:

• Download figure
• Open in new tab
• Download PowerPoint

The objective is to generate a sequence of estimates x^(k) such that the following holds:

$\underset{k\to \infty }{\lim }{x}^{\left(k\right)}=\alpha$ 9

where α ∈ ℝ_n and ∇Φ(α) = 0 (Quarteroni et al., 2010). The sequence {x^(k)} is said to converge to α if, for each pair of successive iterations (x^(k), x^(k+1)), the approximation error decreases according to $\left|x^{(k+1)}-\alpha \right|\le C\left|x^{\left(k\right)}-\alpha \right|$, where C ≥ 0. When the difference between two successive estimates falls below a predefined tolerance or threshold, є, the algorithm is considered to have reached convergence, and thus, the iterations are halted (Quarteroni et al., 2014). Mathematically, this condition can be expressed as follows:

$\left|x^{(k+1)}-x^{\left(k\right)}\right|\le є$ 10

If the threshold is properly chosen, further iterations would not bring significant improvements. The number of iterations required to reach convergence and the actual success of convergence depend on the properties of the objective function, the choice of the initial point x⁽⁰⁾, and the satisfaction of certain regularity conditions (Bierlaire, 2021; Hager, 1988; Nocedal & Wright, 1999; Quarteroni et al., 2014; Yan et al., 2013). These conditions include ensuring that the Jacobian matrix has full rank, that the neglected second-order terms in Equation (7) are either zero or small in magnitude, and that the residuals are small, particularly close to the minimum of the objective function. If any of these conditions are not met, convergence is no longer guaranteed, and the GN method may stagnate, diverge, or converge to a non-stationary point. One commonly employed technique for enhancing the likelihood of convergence is to utilize the GN method in conjunction with a line-search approach. This approach, known as the DGN method (Björck, 1996; Dennis Jr & Schnabel, 1996), involves introducing a multiplicative factor that dampens the Newton step (δx). The correction in Equation (8) is thus replaced by the following:

$x^{(k+1)}=x^{\left(k\right)}+{\mu }^{\left(k\right)}+\delta x^{\left(k\right)}$ 11

The damping factor μ determines the distance traversed along the direction δx. Typically, the value of this factor ranges from 0 to 1. Intuitively, smaller values of μ increase the number of iterations needed for the algorithm to converge, but increase stability, thereby avoiding overshooting the minimum we seek to reach. Conversely, larger values of μ can speed up convergence but may introduce instability that can lead to divergence. When μ = 1, the method reverts to the GN approach in Equation (8).

Theoretically, an optimal μ^(k) value can be found by an exact line-search technique, which involves solving the following minimization problem:

${\mu }^{\left(k\right)}=\underset{\mu \in {ℝ}_0^{+}}{\arg \hspace{0.33em}\min }\Phi (x^{\left(k\right)}+\mu \delta x^{\left(k\right)})$ 12

In the presence of multiple minima, it is common practice to select the first minimum; that is, μ^(k) is the smallest element of $\underset{\mu \in {ℝ}_0^{+}}{\arg \hspace{0.33em}\min }\Phi (x^{\left(k\right)}+\mu \delta x^{\left(k\right)})$ (Bierlaire, 2021).

In practical applications, however, it can be cumbersome to invest extensive effort in solving an additional minimization problem. Instead, an inexact line search based on a trial-and-error approach is often preferred, where different damping factors are tested and the first suitable factor is selected. In the following sections, we will introduce algorithms that implement different methods for selecting the damping factor. Our goal is to demonstrate their potential to enable convergence to the global minimum of the objective function when applied to the position, velocity, and time algorithm based on the LS method. This assessment is particularly important under challenging convergence conditions, such as those encountered in LEO scenarios.

3.1 Simulation Framework

To assess the convergence of the LS method using pseudorange or Doppler shift measurements, we use a method of investigation called “orbit scale down,” which we introduced in a previous work (Morichi et al., 2024). This theoretical method allows the algorithm to be evaluated by simulating variations in the orbital radius of the satellites used for the correction. In the simulations, the satellites were progressively lowered along the line-of-sight (LOS) vector connecting the receiver to the satellite, ensuring that the direction of the receiver-satellite steering vector remains fixed, as shown in Figure 1.

• Download figure
• Open in new tab
• Download PowerPoint

Figure 1

Pictorial view of the orbit scale-down methodology (Morichi et al., 2024); RX: receiver

A total of 100 orbital altitudes were considered in the simulations, with 60 distributed in the LEO region between 160 and 2000 km and 40 in the MEO region up to the orbital altitude of the GPS satellites (20200 km). This approach allowed a trend to be outlined for convergence in the two regions. With each change in the orbital radius of the satellites, both their position and velocity are consistently adjusted. Specifically, by assuming circular orbits, the satellites’ velocity can be adapted based on the orbital period, derived from Kepler’s third law, which is given by the following:

$T_s\approx 2\pi \sqrt{\frac{R_j^3}{GM}}$ 13

where R_j is the orbital radius of the satellite, M is the mass of the Earth, and G is the universal gravitational constant. From Equation (13), the satellite’s velocity for a given orbital radius can be derived. Furthermore, knowing both the initial orbital radius and velocity of the satellite allows us to compute the satellite velocity at a new orbital radius as follows:

$\{{}_{v_{s_2}=\sqrt{\frac{GM}{R_{s_2}}}}^{v_{s_1}=\sqrt{\frac{GM}{R_{s_{{}_1}}}}}\Rightarrow v_{s_2}=v_{s_1}=\sqrt{\frac{R_{s_1}}{R_{s_2}}}$ 14

Given that we are working with synthetic satellite positions and velocities, the measurements used by the LS algorithm must also be adapted. Specifically, both the pseudorange and Doppler shift measurements, originally derived from real observables, must be adjusted accordingly.

To achieve these adjustments, synthetic, ideal measurements were constructed based on the receiver’s ground truth. To make these measurements more realistic, zero-mean Gaussian noise, independent for each satellite, was added. The variance of this noise corresponds to that of the original measurements obtained when the satellites were at their true positions. In the subsequent analyses, we did not model any other types of errors. The focus of this study is to develop an algorithm capable of guiding the estimation of the receiver state toward the global minimum of the objective function, even in scenarios where the traditional GN-based positioning algorithm fails. Provided that the regularity conditions of the GN method are satisfied, the presence of errors and distortions in the measurements may affect the rate of convergence but does not hinder overall convergence. The minimum is still reached, although it may be shifted from its nominal location. As a result, omitting these errors does not compromise the generality of our findings, as any convergence issues observed can be attributed solely to the algorithm and the geometry of the analyzed scenario.

The effectiveness of the proposed algorithm was initially tested using a static receiver located in Turin, Italy (45.0653˚ N, 7.6589˚ E, 321.11 m). We started with real data collected through an antenna and front-end (i.e., Ettus Research USRP N210), which were then processed with a software receiver implemented in MATLAB. The real scenario involved the use of five GPS satellites for the correction. Based on the position and velocity information of these five satellites, obtained via post-processing of the collected data, we applied the orbit scale-down framework presented above. This approach allowed us to calculate the position and velocity of the same five satellites at different orbital altitudes. Unlike the approach described by Morichi et al. (2024), where the algorithm’s convergence relied on gradually adjusting the initialization points of the GN method as the satellite’s orbital altitude decreased, in this case, no prior knowledge about the receiver position is assumed. Consequently, the initial estimate of the receiver state is set to all zeros, as in conventional approaches, meaning that the position state is placed at the center of the Earth and all other components related to the bias (or drift) of the receiver clock are assumed to be zero.

3.2 ADGN Algorithm

To enhance the resilience of the GN-based positioning algorithm to convergence issues, we propose utilizing the DGN method. This section presents the algorithm that we have developed, which is based on the theoretical framework outlined in Section 2. The pseudocode for the ADGN method, which is our adaptation of the GN method, is shown in Algorithm 1. This method is adaptive in that it introduces a damping factor, μ, that is dynamically adjusted: the damping factor increases when the algorithm converges and decreases when the algorithm does not converge.

To determine whether the estimates are converging between iterations, it is not possible to use the error estimate as presented in Section 2.2, because the true value of the minimum is not a priori known. However, taking advantage of the fact that we are looking for the minimum of the objective function, it is possible to check the descent condition between two successive iterations. Specifically, the goal is to find a value for the damping factor μ that satisfies the following:

$\Phi \underset{x^{(k+1)}}{\underbrace{(x^{\left(k\right)}+{\mu }^{\left(k\right)}\delta x^{\left(k\right)})}}\le \Phi \left(x^{\left(k\right)}\right)\Rightarrow \underset{x^{(k+1)}}{\underbrace{\Vert R(x^{\left(k\right)}+{\mu }^{\left(k\right)}\delta x^{\left(k\right)})\Vert }}\le \Vert R\left(x^{\left(k\right)}\right)\Vert$ 15

This formula indicates whether the objective function, computed using the updated state estimate, x^(k+1), decreases. The demand for exact equality may be overly stringent, leading to unnecessary iterations when the difference between two consecutive iterations of the objective functions is negligible. Specifically, we calculate the tolerance as $\left|R\left(x^{(k+1)}\right)-R\left(x^{\left(k\right)}\right)\right|<\tau$, where τ is an arbitrary level of tolerance allowing the descent condition to be validated even when the residuals of two successive iterations are not perfectly equal. If Equation (15) is satisfied, the algorithm is considered to be converging, and in this case, the damping factor μ is increased for the next iteration to accelerate the descent to the minimum. Conversely, if Equation (15) is not met, an inexact line-search method is employed. The damping factor μ is iteratively halved until a value that satisfies the condition in Equation (15) is found. The results presented in the following sections employ a tolerance of τ = 10^–8.

The algorithm continues iterating until convergence is achieved. Convergence is determined by monitoring the norm of the applied Newton step, $\Vert {\mu }^{\left(k\right)}\delta x^{\left(k\right)}\Vert$, which indicates the magnitude of the correction applied to the current state estimate. Once this norm becomes sufficiently small, indicating that further iterations will produce negligible changes, the algorithm terminates. In our simulations, the threshold for this check was set using the parameter є = 10^–5 for both the pseudorange- and Doppler-based algorithms.

To mitigate problems associated with the two iterative loops of the algorithm, it is necessary to introduce two exit conditions. The first condition halts the convergence iteration loop if the number of iterations exceeds a predefined limit, set to $\zeta =500$ in our simulations. The second condition concludes the routine for determining the appropriate damping factor, μ, and is triggered when μ falls below a specified threshold, η, which is set to 10^–6 in our simulations. If the algorithm exits based on either of these conditions, the simulation is considered non-converging. This outcome suggests that no value of μ leads to a decrease in the objective function, indicating that the algorithm may be following an incorrect direction.

3.2.1 ADGN Algorithm Results

The ADGN algorithm was tested with both pseudorange and Doppler shift measurements, with the same experiments performed separately for each measurement type. The theoretical “orbit scale-down” framework was applied to both positioning algorithms, enabling an evaluation of the performance with satellites at various orbital altitudes. For each orbital radius, the algorithm was initialized with eight different damping factors ( μ ∈{0.01,0.02,0.04,0.08,0.16,0.32,0.64,1}), and the average number of iterations required for convergence across epochs was monitored. Additionally, we tested a different adaptive method for the damping factor, referred to as 1B in the following figures, which employs a backtracking strategy by setting the damping factor to 1 in line 6 of Algorithm 1. This approach ensures that the damping factor does not retain memory of the previous iteration’s state; the damping factor always starts from 1 and only decreases when Equation (15) is not met. To ensure that the algorithm converged on the desired minimum, the position found was compared with the receiver position ground truth.

Figure 2 presents the baseline reference results obtained from the LS method solved by the GN method. The heatmap shows the orbital altitude of the satellites on the x-axis and versions of the algorithm based on pseudorange and Doppler shift on the y-axis. The color gradient represents the average number of iterations across the epochs for each scenario. The MEO and LEO regions are graphically marked by a vertical split on the heatmaps. The black region highlights scenarios in which the algorithm failed to converge. This failure occurred predominantly in the LEO region. The presented results establish a benchmark for interpreting and gaining a better understanding of the advantages of employing the proposed technique.

• Download figure
• Open in new tab
• Download PowerPoint

Figure 2

Average number of iterations for the LS method solved via the GN method

Figure 3 presents the results of the ADGN algorithm. Each heatmap depicts satellite orbital altitudes on the x-axis, while the y-axis displays various initializations of the damping factor μ, with colors indicating the average iteration count. Figures 3(a) and 3(b) show the convergence results for the ADGN algorithm, applied to pseudorange and Doppler shift measurements, respectively. As in Figure 2, scenarios in which the algorithm fails to converge are shown in black.

• Download figure
• Open in new tab
• Download PowerPoint

Figure 3

Average number of iterations for the LS method solved using the ADGN algorithm presented in Algorithm 1 (a) Convergence heatmap (pseudorange) (b) Convergence heatmap (Doppler shift)

Figures 3(a) and 3(b) both show a clear trend. Low values for the initial damping factor μ⁽⁰⁾ ensure convergence for each satellite altitude considered. This result arises because the correct descent path can be followed from the first iteration, with the damping factor progressively doubled until convergence is achieved. A higher damping factor in the initial iterations does not always guarantee convergence, as it can cause the algorithm to overshoot the minimum from the very start. Additionally, the heatmaps in Figure 3 indicate that higher damping factor values lead to faster convergence for satellites in higher orbits, while the number of iterations required increases as satellite altitude decreases. This trend is observed because, to satisfy Equation (15), the damping factor must first decrease, necessitating more iterations to return to 1. To address this issue, we also tested an alternative adaptive method, where the damping factor is reset to 1 at each new iteration. However, although this approach records fewer iterations in some successful cases, returning to 1 at every iteration can introduce some risk, resulting in a significantly higher number of divergence cases compared with other initializations. By averaging the number of iterations needed for convergence across all satellite altitudes, we observed that the damping factor yielding the fewest average iterations across all satellite altitudes is 0.08 for pseudorange, with an average iteration count of 9.17. Conversely, for the Doppler shift, a damping factor value of 0.02 is found, with an average iteration count of 11.5.

• Download figure
• Open in new tab
• Download PowerPoint

ALGORITHM 1

ADGN Algorithm

The upcoming subsections propose modifications to the ADGN algorithm to enhance its performance and guarantee convergence for all orbital radii and damping factor initializations.

3.3 AW-DGN Algorithm

The descent condition presented in Equation (15) is not always sufficient. In fact, it is well known in the mathematical literature that the mere verification of Equation (15) is not sufficient to confirm the appropriateness of the damping factor (Bierlaire, 2021; Quarteroni et al., 2014). When choosing an adequate damping factor, it is important to remember that the damping factor cannot be too large or too close to zero (Bierlaire, 2021). Whereas higher values of the damping factor may reduce the number of iterations required for convergence, excessively large values risk overshooting the desired minimum, resulting in non-convergence or convergence to a different minimum. Conversely, a damping factor too close to zero can completely halt progress, which is equally detrimental.

Therefore, criteria that are stricter than Equation (15) are often employed, such as Wolfe conditions (Dennis Jr & Schnabel, 1996; Gill et al., 2019). These conditions establish an upper and lower bound for the dimension of the damping factor, guiding the selection of an appropriate size, and can be expressed as follows:

$\Phi (x^{\left(k\right)}+\mu \delta x^{\left(k\right)})<\Phi \left(x^{\left(k\right)}\right)+\mu {\beta }_1\nabla \Phi {\left(x^{\left(k\right)}\right)}^{\top }\delta x^{\left(k\right)}$16a

$\nabla \Phi {(x^{\left(k\right)}+\mu \delta x^{\left(k\right)})}^T\delta x^{\left(k\right)}\ge {\beta }_2\nabla \Phi {\left(x^{\left(k\right)}\right)}^T\delta x^{\left(k\right)}$16b

where β₁ and β₂ are two parameters, usually chosen such that 0 < β₁ < β₂ < 1. Once these parameters are determined, the objective is to find a damping factor μ that satisfies both conditions.

The first Wolfe condition (Equation (16a)), also known as Armijo’s rule (Dennis Jr & Schnabel, 1996; Quarteroni et al., 2014), is based on the concept of a sufficient decrease in the objective function between two successive iterations. The use of this condition aims to prevent the selection of damping factors that are large but fail to result in a significant decrease in the objective function. The decrease in objective function should be proportional to the slope of the function in x^(k) along the direction δx^(k).

In contrast, the second Wolfe condition (Equation (16b)), also known as the curvature condition, prevents the damping factor from being too small by checking the slope at x^(k) + μδx^(k) and at x^(k).

Figure 4 provides a two-dimensional pictorial representation, inspired by Quarteroni et al. (2014), that helps clarify the meaning of the Wolfe conditions. For the purpose of this example, we consider a univariate objective function $\Phi \left(\mu \right)=\Phi (x^{\left(k\right)}+\mu \delta x^{\left(k\right)})$. The goal is to find a minimum of this function, ideally the global minimum. In Figure 4, Φ(μ) is shown as a solid black curve. At iteration k, the function can be approximated using a first-order Taylor expansion around x^(k), i.e., $\Phi \left(x^{\left(k\right)}\right)+\mu \nabla \Phi \left(x^{\left(k\right)}\right)\delta x^{\left(k\right)})$, which is shown as a dashed blue line in Figure 4(a). The first Wolfe condition requires that, after a step of size μ has been taken, the objective function value Φ(μ) lie below this linear approximation. The values of μ that satisfy this condition are highlighted in blue in Figure 4(a). The second Wolfe condition concerns the slope of Φ(μ) Specifically, this condition ensures that the derivative at the chosen μ is not smaller than the slope at x^(k), thereby preventing the selection of excessively small damping factors. In Figure 4 b, the dashed blue line represents ∇Φ(x^(k)) δx^(k) and the shaded regions indicate the values of μ for which the second Wolfe condition is satisfied. Finally, in Figure 4(c), the points along the curve Φ(μ) highlighted in blue represent the values of μ that simultaneously satisfy both Wolfe conditions.

• Download figure
• Open in new tab
• Download PowerPoint

Figure 4

Pictorial representation of the Wolfe conditions for a univariate objective function Φ(μ) (a) First Wolfe condition (b) Second Wolfe condition (c) Combined view of both Wolfe conditions

The points highlighted in blue correspond to points that satisfy the Wolfe conditions.

In scenarios with noise-affected measurements, such as the scenario under consideration, applying Wolfe conditions to guide the choice of a damping factor can be challenging. Additionally, these conditions depend on two parameters, β₁ and β₂, which must be set according to the specific scenario. While the mathematical literature suggests values for these parameters (Bierlaire, 2021; Quarteroni et al., 2014) (typically β10^–4 and β₂ = 0.99), experiments conducted by the authors have shown that these values are often ineffective. Generalizing parameters that work across various scenarios is not straightforward. Although an adaptive approach for adjusting β₁ and β₂ during iterations may be preferable, it too often proves to be ineffective. Furthermore, computing Equation (16b) requires information about the prime derivatives at the point UT, adding computational burden to the problem. For these reasons, we aim to leverage the principles behind Wolfe conditions and propose a heuristic solution that is more computationally efficient in terms of parameter settings and cost.

We first focus on the first Wolfe condition, which aims to ensure a sufficient decrement in the objective function, proportional to the damping factor. However, because a numerical method is employed, the specific context of the problem to which the condition applies is not considered, but the context is relevant and can be exploited. In our scenario of interest, we are dealing with a positioning problem, where the primary objective is to obtain an accurate estimate of the receiver states.

The knowledge that the receiver is situated on Earth can be exploited to our advantage and forms the basis of our heuristic approach. Our idea is to prevent the iterative approximation point from deviating too far from the Earth’s surface. This strategy helps us avoid local minimum points located off the Earth’s surface, which are unlikely to correspond to the receiver’s actual position. Therefore, in addition to the descent condition in Equation (15), we introduce an additional criterion that triggers a decrease in damping factor even when the estimated receiver position extends beyond the Earth’s surface.

In terms of pseudocode, the modifications are minimal and involve adjusting the verification of the descent condition (line 9 of Algorithm 1) as follows:

$\Vert R\left(x^{UT}\right)\Vert \Vert <\Vert R(x)\Vert \wedge \Vert {\{x\}}_{i=1}^3\Vert <\psi +\theta$ 17

where ${\{x\}}_{i=1}^3$ represents the components of the receiver state x related to its position, ψ identifies the mean radius of the Earth (ψ = 6.371 km), and θ is a tolerance parameter. In our simulations, θ has been set to 50 km, but this term can be tailored to specific applications. The ADGN algorithm with this condition added will be referred to as AW-DGN from here on.

This condition is beneficial not only for implementation purposes, as it does not require any additional calculations beyond what is already available, but also for its operational efficiency. This condition is particularly relevant during the initial iterations, when there is a risk of the position estimate drifting too far away, which could cause the algorithm to diverge or converge to other local minima. After properly adjusting the initial steps, the algorithm progresses smoothly towards the minimum point, and the second condition no longer applies. This approach ensures that the algorithm is not overly restrictive and allows it to converge effectively.

The second Wolfe condition (Equation (16b)), which prevents choosing a damping factor too close to zero, is not explicitly verified. This is because the proposed algorithm already imposes a constraint on the damping factor, ensuring that it cannot be smaller than η, which is set to 10^–6 in our simulations. Considering that the initialization of the damping factor used is at least four orders of magnitude larger than η and that the algorithm gradually reduces the damping factor whenever the descent condition is not satisfied, the chosen damping factor is inherently chosen not to be too small, making our implementation significantly simpler compared with enforcing the second Wolfe condition (Equation (16b)).

3.3.1 AW-DGN Algorithm Results

Figure 5 shows two heatmaps for the LS method based on pseudorange and Doppler shift measurements. The color scale indicates the average number of iterations across epochs required by the algorithm to reach convergence. For each satellite orbital altitude, nine different damping factor initializations were tested, including the backtracking approach labeled as “1B” in the heatmap. A notable difference between Figure 5 and Figure 3 is the absence of black areas that indicate convergence failure in the algorithm. This result demonstrates that the modifications applied to the algorithm have been effective, resulting in convergence for every scenario analyzed.

• Download figure
• Open in new tab
• Download PowerPoint

ALGORITHM 2

AW-DGN-BT

• Download figure
• Open in new tab
• Download PowerPoint

Figure 5

Average number of iterations for the LS method solved using the AW-DGN algorithm (a) Convergence heatmap (pseudorange) (b) Convergence heatmap (Doppler shift)

Furthermore, upon closer examination, it is evident that for a lower damping factor, the number of iterations remains unchanged compared with the previous simulation, as the applied modification does not directly affect the number of iterations. However, for higher damping factors, where the modification of the descent condition comes into play, differences can be observed. Nevertheless, the general trend observed previously is confirmed, wherein for higher initialization damping factor values, the number of iterations increases as the orbital altitude of the satellites decreases.

The damping factor that resulted in the fewest average iterations across all sets of orbital radii was the 1B initialization for both the pseudorange and Doppler shift, with average iteration counts of 7.23 and 7.5, respectively. This result suggests that by using Equation (17), which prevents algorithm divergence, the damping factor with backtracking to 1 at each iteration offers the best performance in terms of convergence iterations compared with all other damping factor initializations, especially in the LEO region. Thus, this approach is the most effective among all tested implementations. Algorithm 2 presents the pseudocode with all modifications needed to implement the AW-DGN with backtracking (AW-DGN-BT).

4.1 Simulation Framework

We simulated a Walker Delta constellation similar to the Galileo system, consisting of 24 MEO satellites. Using MATLAB’s built-in orbit propagator, we computed the satellite positions and velocities over a 10-h period, with sampling every 10 s, which resulted in 600 epochs. The receiver position was simulated at three locations: Cape Town, Longyearbyen, and Shanghai. A visibility analysis was then conducted to determine which satellites were usable for each location. The large time window provides a significant advantage in terms of satellite geometry, as the satellite positions change not only between different scenarios but also within the same scenario over time. This framework incorporates the orbit scale-down approach described in Section 3.1. For each simulated location, both satellite positions and velocities were adjusted, and the positioning algorithm (Algorithm 2) was tested to evaluate its robustness in terms of convergence. Table 1 summarizes the three scenarios analyzed.

4.2 Simulation Results

Figure 6 compares the average number of iterations across epochs required for convergence among the GN-based positioning algorithm, ADGN-BT, and AW-DGN-BT, across the different scenarios summarized in Table 1. Each column in the heatmaps corresponds to a different satellite altitude, whereas rows distinguish between the three algorithms being compared. For each scenario, both pseudorange and Doppler measurements were tested. Black areas in the heatmaps indicate cases where the algorithm did not converge, which occur frequently for both the GN and ADGN-BT methods as the satellite orbital radius approaches the LEO region. In contrast, the AW-DGN-BT algorithm is highly effective, consistently achieving convergence in all scenarios. These trends, visually illustrated in Figure 6, are quantitatively summarized in Tables 2 and 3. In particular, Table 2 reports the percentage of successful convergence for each algorithm and measurement type across the three representative scenarios. The AW-DGN-BT algorithm achieves a consistent 100% success rate in all cases, unlike GN and ADGN-BT, which show significantly lower success percentages (below 40% in virtually all scenarios). Table 3 further reports the average number of iterations and execution time (in milliseconds) across all considered satellite altitudes and scenarios (Table 1) for each algorithm and measurement type. These averages are calculated only over the satellite altitudes where the algorithm successfully converged. Because AW-DGN-BT converges in all cases, its averages are computed over all 100 satellite altitudes, whereas those for GN and ADGN-BT include only the subset of altitudes that resulted in successful convergence. Although AW-DGN-BT requires slightly more iterations (around 7) and marginally higher computation time (around 1 ms), this trade-off is offset by its superior convergence reliability. Note that the execution time is reported here as a representative metric, although it may vary depending on implementation details, hardware specifications, and the programming environment used.

• Download figure
• Open in new tab
• Download PowerPoint

Figure 6

Average number of iterations for the LS method solved using both the GN and AW-DGN-BT algorithms; LLA: latitude, longitude, altitude (a) S1a: Cape Town (South Africa), LLA = [–33.92899°, 18.41740°, 0 m] (b) S1b: Longyearbyen (Norway), LLA = [78.22894°, 15.39493°, 0 m] (c) S1c: Shanghai (China), LLA = [31.14274°, 121.80411°, 5 m]