George Mason University

Project 2: GPS, Conditioning, and Nonlinear Least Squares

The goal of this project is to calculate the position of a global positioning system (GPS) receiver using the information broadcast by the GPS satellites. The 24 satellites of the GPS system carry atomic clocks to relay a time stamp as well as their position in space. This transmission for a satellite \( i \) contains the four parameters \( (x_i, y_i, z_i, t_i) \).

When the time stamp arrives at the receiver it is easily converted into the transmission travel time by subtracting the current time on the receiver from the satellite's time stamp. This is then converted into the satellite's distance from the receiver by multiplying the travel time by the speed of light, \( c = 299792.458 \frac{\text{km}}{\text{s}} \).

However, while the atomic clocks on board the satellites are accurate to \( 10^{-9} \) seconds per day and are synchronized regularly, the clock of a typical GPS receiver does not have this precision. This error in time synchronization, denoted \( d \) must also be solved for if an accurate position is to be computed using the GPS.

As the goal will be to solve for the position as \( (x, y, z, d) \), for clarity's sake each satellite \( i \) will contain the same parameters labelled as \( (A_i, B_i, C_i, t_i) \). The solution to the position relies on the fact that three satellites only intersect at two points.

One of these points will be near the earth and relevant to the purposes of the position problem. With the inclusion of the time error, a system of four equations using four satellites can be used to solve for the position and time error. This is described by the system directly below.

System 1.1

\[ r_1(x, y, z, d) = \sqrt{ (x - A_1 )^2 + (y - B_1)^2 + (z - C_1)^2 } - c(t_1 - d) = 0 \\ r_2(x, y, z, d) = \sqrt{ (x - A_2 )^2 + (y - B_2)^2 + (z - C_2)^2 } - c(t_2 - d) = 0 \\ r_3(x, y, z, d) = \sqrt{ (x - A_3 )^2 + (y - B_3)^2 + (z - C_3)^2 } - c(t_3 - d) = 0 \\ r_4(x, y, z, d) = \sqrt{ (x - A_4 )^2 + (y - B_4)^2 + (z - C_4)^2 } - c(t_4 - d) = 0 \]

As five to eight satellites are in direct line of sight from any point on earth, the inputs of this system should be available to a working GPS receiver.

[Return to Top]

Step 1:

We employed the Multivariate Newton's Method to solve for the position and time correction using four satellites.

The position solution for \((x, y, z, d)\).

We will test the Multivariate Newton's Method for 4 satellites \((A_i, B_i, C_i, t_i)\) with a known answer to see if we can obtain the same solution.

System 1.1 is algebraically manipulated into the form of System 1.2 below.

\begin{align*} f_1 = (x - A_1)^2 + (y - B_1)^2 + (z - C_1)^2 - (c(t_1 - d))^2 = 0\\ f_2 = (x - A_2)^2 + (y - B_2)^2 + (z - C_2)^2 - (c(t_2 - d))^2 = 0\\ f_3 = (x - A_3)^2 + (y - B_3)^2 + (z - C_3)^2 - (c(t_3 - d))^2 = 0\\ f_4 = (x - A_4)^2 + (y - B_4)^2 + (z - C_4)^2 - (c(t_4 - d))^2 = 0 \end{align*}

To solve for \(\alpha = (x, y, z ,d)\) the the Multivariate Newton Method takes the Jacobian matrix of System 1.2 as \(F(\alpha) = (f_1, f_2, f_3, f_4)\) as seen below.

\[\text{DF}(\alpha) = \left( \begin{array}{cccc} \frac{\partial f_1}{\partial x} & \frac{\partial f_1}{\partial y} & \frac{\partial f_1}{\partial z} & \frac{\partial f_1}{\partial t} \\ \frac{\partial f_2}{\partial x} & \frac{\partial f_2}{\partial y} & \frac{\partial f_2}{\partial z} & \frac{\partial f_2}{\partial t} \\ \frac{\partial f_3}{\partial x} & \frac{\partial f_3}{\partial y} & \frac{\partial f_3}{\partial z} & \frac{\partial f_3}{\partial t} \\ \frac{\partial f_4}{\partial x} & \frac{\partial f_4}{\partial y} & \frac{\partial f_4}{\partial z} & \frac{\partial f_4}{\partial t} \end{array} \right) = \left( \begin{array}{cccc} 2(x - A_1) & 2(y - B_1) & 2(z - C_1) & 2c^2(t_1 - d) \\ 2(x - A_2) & 2(y - B_2) & 2(z - C_2) & 2c^2(t_2 - d) \\ 2(x - A_3) & 2(y - B_3) & 2(z - C_3) & 2c^2(t_3 - d) \\ 2(x - A_4) & 2(y - B_4) & 2(z - C_4) & 2c^2(t_4 - d) \end{array} \right)\]

Beginning with an initial vector \(\alpha_0\) we iterate through \(n \in \Bbb N\) steps of the process of solving

\[\left\{ \begin{array}{lr} DF(\alpha_k)s = -F(\alpha_k) \\ \alpha_{k+1} = \alpha_k + s \end{array} \right. \hspace{1 cm} \text{for } k = 0, 1,\ldots , n\]

Depending on the initial vector, this iteration will converge to one of two possible solutions.
To solve for the far from Earth solution, we used the initial input \((x, y, z, d) = (1000, 1000, 50000, 0) \).
Likewise, to solve for the solution on the Earth, we used the initial input \((x, y, z, d) = (0, 0, 6370, 0)\).
Figure 1.1 below displays the two solutions achieved using these initial vectors.

[Return to Top]

Step 2:

In order to solve System 1.2 as a closed form solution each component \(x\),\(y\), and \(z\) were isolated in terms of \(d\). By subtracting the last three equations from the first, three linear equations can be obtained as shown below in System 2.1.

System 2.1:

\begin{align*} 2(A_4 - A_1)x + 2(B_4 - B_1)y + 2(C_4 - C_1)z - 2c^2(t_4 - t_1)d - c^2(t_1^2 - t_4^2) + A_1^2 - A_4^2 + B_1^2 - B_4^2 + C_1^2 - C_4^2 = 0 \\ 2(A_3 - A_1)x + 2(B_3 - B_1)y + 2(C_3 - C_1)z - 2c^2(t_3 - t_1)d - c^2(t_1^2 - t_3^2) + A_1^2 - A_3^2 + B_1^2 - B_3^2 + C_1^2 - C_3^2 = 0 \\ 2(A_2 - A_1)x + 2(B_2 - B_1)y + 2(C_2 - C_1)z - 2c^2(t_2 - t_1)d - c^2(t_1^2 - t_2^2) + A_1^2 - A_2^2 + B_1^2 - B_2^2 + C_1^2 - C_2^2 = 0 \end{align*}

Note that System 2.1 can be broken into column vectors in terms of the variables \(x, y, z, d,\) and constants by first taking the column vectors of System 2.1 and denoting them as \(x \vec{u}_x, y \vec{u}_y, z \vec{u}_z, d \vec{u}_d,\) and a constant term \(\vec{w}\). This can be rewritten as

\(x \vec{u}_x + y \vec{u}_y + z \vec{u}_z + d \vec{u}_d + \vec{w} = \vec{0}\).

A formula for \(x\) in terms of \(d\) can be acquired from

\(\vec{0} = det[\vec{u}_y | \vec{u}_z | x \vec{u}_x + y \vec{u}_y + z \vec{u}_z + d \vec{u}_d + \vec{w}] \)

Likewise, a formula for \(y\) and \(z\) in terms of \(d\) can be acquired from

\(\vec{0} = det[\vec{u}_x | \vec{u}_z | x \vec{u}_x + y \vec{u}_y + z \vec{u}_z + d \vec{u}_d + \vec{w}] \)
\(\vec{0} = det[\vec{u}_x | \vec{u}_y | x \vec{u}_x + y \vec{u}_y + z \vec{u}_z + d \vec{u}_d + \vec{w}] \)

Finally, by substituting \(x, y, \text{and } z\) in terms of \(d\) back into row one of System 1.2, we obtain a quadratic equation in terms of \(d.\) The quadratic equation, popularly represented as below, is then used to solve for both roots of \(d.\)

\[r = \frac{-b \: \pm \: \sqrt{b^2 - 4ac}}{2a} \]

The results of this method are shown in Figure 2.1, and the associated code can be viewed by clicking on the table.

Figure 2.2 shows the difference in position and time synchronization error for the Multivariate Newton Method and the closed form solution. There was little difference between the two; the scale of the differences in position were in micrometers.

[Return to Top]

Step 4:

We generalized our Multivariate Newton's method procedure from Part 1 into a function capable of estimating the position of an observer on the north pole based on the signal from any set of 4 satellites.

Unlike previous examples, which specified satellite positions as a set of \((x, y, z)\) coordinates, our function for this part took as its arguments spherical coordinates \(\phi\) (ranged 0 to \(\pi\)) and \(\theta\) (ranged 0 to \(2\pi\)). All satellites were assumed to be 20km above the earth's surface \((\rho =26570 \text{km})\). As shown below, we had MATLAB generate a 3D plot to roughly show the positioning of the satellites.

In order to determine the sensitivity of the measured coordinates to errors in the timing of the input data, the algorithm applied tiny changes (\(\pm d\) for some small \(d\), such as 10 nanoseconds) to the observed time delays and measured the largest coordinate error in the computed result.

Figure 4.1. Click on the image to view the associated code.

All possible combinations of deviations (\(-d\) error all four time delays, \(+d\) error for all four time delays, no error, any combination thereof), a total of \(3^4=81\) possibilities, were tried in order to determine the maximum output error, and the unit-less error magnification factor (EMF) was calculated as follows:

\[\text{emf} = \frac{|| \Delta x, \Delta y, \Delta z ||_{\infty}}{c || \Delta t_1 \Delta t_2 \Delta t_3 \Delta t_4 ||_{\infty}} = \frac{|| \Delta x, \Delta y, \Delta z ||_{\infty}}{cd}\]

Example:
>> part4(10^-8, [0, 0.4, 0.8, 1.2], [0, pi/2, pi, 3*pi/2])
maxerr =
 0.0132
ans =
 4.3946

For the given coordinates, we have a maximum error magnification factor of 4.3946, meaning that a timing error of 10ns (2.998 meters backwards error, the distance that the GPS signal travels in that time) translates to a forward error of 13.2 meters. Our estimate for the condition number of the problem at these coordinates is 4.3946.

[Return to Top]

Step 5:

In Step 4 we chose widely spaced satellite coordinates, and got a condition number of 4.3946. We will now proceed to alter the spacing of the satellites in order to check the error and condition number under differing conditions.

Figure 5.1. Click on the image to view the associated code.

Case 1: When the satellites are somewhat clustered.

This is done first without error, and then with 10 nanoseconds of error:

>> part4(0, [0.5, 0.75, 1.0, 1.25], [1.5, 2.0, 2.5, 3.0])
maxerr =
 3.2742e-11
ans =
 Inf

>> part4(10^-8, [0.5, 0.75, 1.0, 1.25], [1.5, 2.0, 2.5, 3.0])
maxerr =
 0.2149
ans =
 71.6852

Below is the generated 3D plot of Case 1, Figure 5.1.

As shown by the first function call (with no input error), the equation is still solvable at machine precision (with an output error in the nanometers), the error magnification factor is much higher than before; in this instance, a timing error of 10 nanoseconds resulted in an error of 215 meters for the calculated position.

Case 2: When the satellites are tightly clustered.

Figure 5.2.

Once again, the code is run first without error, then with 10 nanoseconds of error:

>> part4(0, [0.70, 0.71, 0.72, 0.73], [2.15, 2.1, 2.05, 2.0])
maxerr =
 1.0927e-07
ans =
 Inf

>> part4(10^-8, [0.70, 0.71, 0.72, 0.73], [2.15, 2.1, 2.05, 2.0])
maxerr =
 859.5457
ans =
 2.8671e+05

Below is the generated 3D plot of Case 2, Figure 5.2.

A Special Case: When all the satellites are at an identical inclination

Figure 5.3.

When \(\phi\) is the same for all four satellites, the method reports a singular matrix and fails completely, despite the satellites being spaced equally far apart:

>> part4(10^-8, [0.5, 0.5, 0.5, 0.5], [0, pi/2, pi, 1.5*pi])

Below is the generated 3D plot of this case, Figure 5.3.

Interestingly enough, this does not occur when the value of \(\theta\) is held constant. We determined that this is because the observer is at the north pole. Thus if all four satellites are at the same latitude, they are equidistant from the observer and depending on the value of d could place the observer at any point along the z-axis.

[Return to Top]

Step 6:

The goal of step six is to discover whether the GPS error and condition number can be reduced by increasing the number of satellites used to compute the position and time synchronization error. As there are four variables being solved for, \(x, y, z, \text{ and }d\), the addition of one or more satellites will not allow our matrix to be square any longer. We employ the Gauss-Newton Method to handle this scenario.

We will briefly describe the process of the Gauss-Newton method. If \(n\)satellites are employed for the solution, with \(n \geq 4\), then a matrix with dimensions \(n \times 4\) is produced.

This matrix,

\( \left( \begin{array}{c} r_1(x, y, z, d)\\ r_2(x, y, z, d)\\ r_3(x, y, z, d)\\ \vdots\\ r_n(x, y, z, d) \end{array} \right)\)

will then have its Jacobian taken as \( A = Dr(x, y, z, d) \).

Denote \(\alpha_0 = ( x_0, y_0, z_0, d_0 )\) We converge to a solution by iterating through \(m \in \Bbb N \) steps of solving

\(\left\{ \begin{array}{lr} A^TAv_k = -A^Tr(\alpha_k) \\ \alpha_{k+1} = \alpha_k + v \end{array} \right. \hspace{1 cm} \text{for } k = 0, 1, \ldots , m\).

We used \(\alpha_0 = (0, 0, 6370, 0)\) as our initial vector. Figure 6.1 below displays the condition number of each matrix \(A^TA\) and the position error calculated as \(\sqrt{x_{err}^2 + y_{err}^2 + z_{err}^2}\).

Figure 6.2. Click on the image to view the associated code.

Figures 6.1 and 6.2 show that increasing satellites reduces the EMF and provides a reduction in error. While the number of satellites required for increased precision is high, there is a benefit associated with them in the Gauss-Newton Method, though this benefit gets smaller and smaller the more you add.

[Return to Top]

Conclusion:

Figure 7.1.

In Steps One and Two, the Multivariate Newton Method and a closed form solution using the quadratic formula were employed to check the validity and accuracy of the solutions they output. Both procedures produced nearly identical results against a provided satellite arrangement and solution with little trouble.

Step Four and Step Five saw the positions selected in order to check the conditioning of the procedures. In this case the Multivariate Newton's Method was employed.

The three examples addressed showed a significant increase in error magnification factor with more clustered satellite arrangements. The well-spaced-out satellites in Step Four had an EMF of 4.3946, while the moderately clustered satellites in Step Five, Case 1 had an EMF of 71.6852, and the nearly-touching satellites in Step 5, Case 2 had an EMF of 28,671.

To try to investigate this further, a script was constructed to generate random sets of satellites and compare the angular distance between them and the EMF of the results they produced. While it is obvious that there are numerous factors involved, one can see that the EMF is roughly proportionate to the angular distance in Figure 7.1.

Finally, Step Six set out to reduce error by adding more satellites. The Gauss Newton method has no issue with clustered satellites being part of the array so long as adequately spaced out satellites are available. As sattelites are added, the EMF is reduced and no loss of precision occurs. The Gauss Newton method ensures that, by using the maximum available satellites in the GPS system, clustering issues can be avoided and error can be reduced significantly.

[Return to Top]