The global positioning system or GPS is a constellation of satellites that are used to approximate the location of a GPS reciever on Earth. GPS uses 4 satellites that are in view of the reciever to solve 4 equations for the \((x,y,z)\) coordinates of the reciever and a value \(d\) which is the difference in time between the reciever's clock and the satellites' clocks. The difference in time is important because the location is approximated using the travel time of a signal from the satellite to the reciever, and this happens in much less than a second. The difference in time can cause error because the satellites' clocks are accurate to \(10^{-8}\) of a second while the reciever's clock is much less accurate.
The equations that are solved to approximate a reciever's location using GPS are: \begin{align*} (x-A_1)^2+(y-B_1)^2+(z-C_1)^2-(c(t_1-d))^2&=0 \\ (x-A_2)^2+(y-B_2)^2+(z-C_2)^2-(c(t_2-d))^2&=0 \\ (x-A_3)^2+(y-B_3)^2+(z-C_3)^2-(c(t_3-d))^2&=0 \\ (x-A_4)^2+(y-B_4)^2+(z-C_4)^2-(c(t_4-d))^2&=0 \end{align*} where \(x, \ y, \text{ and } z \) are the rectangular coordinates of the reciever, \(A,\ B, \text{ and } C \) are the coordinates of the satellites, \(d\) is the difference in time between the reciever and the satellite's clocks, and \(t\) is the travel time for the signal from the satellite to the reciever. We will demonstrate the approximation of a reciever's location by solving the system of equations for specific satellite coordinates using the quadratic equation and Newton's method. We will then demonstrate sources for error in the approximation by choosing different satellite positions and amounts.
Using the equations listed in the introduction, we will approximate the location of a reciever. In order to use the equations to determine the location of the reciever, an equation solving method is needed. In this case, we will use the multivariate Newton's method.
The multivariate Newton's method iterates through successive guesses until one of the guesses is close enough to the correct answer to be a good approximate guess. To find each guess, the multivariate Newton's method solves the linear system \[ Dg=-F \] where \(F\) is the vector of the function at the current guess, \(D\) is the Jacobian matrix of the function at the current guess, and \(g\) is the vector of unknown variables which will be the next guess. This method ususlly converges for the GPS equations to 7 correct digits in less than 20 iterations. Our code for the multivariate Newton's method used to solve the GPS equations can be found in the file MVNewtonsInputs. Our Newton's method program also uses the functions FCreator and DCreator to create the function and its Jacobian at the current guess.
Given the satellite coordinates \( (15600,7540,20140),(18760,2750,18610),(17610,14630,13480),(19170,610,18390) \), the time intervals \( (0.07074,0.07220,0.07690,0.07242) \) and using the initial guess \((0,0,6370,0)\), we use the multivariate Newton's method to approximate the position of the reciever of the GPS signals. In this case the approximate coordinates of the reciever are \[ (x,y,z,d)=(-41.772709, -16.789194, 6370.059559, -0.003201) \] where \(x,y,z\) are the coordinates and \(d\) is the difference in time between the reciever and the satellite's clocks.
There is a second solution to the GPS equations for the same inputs as above, but the second solution is a solution where to coordinates are off of Earth. Using the same inpus as above, the second solution to the GPS equations is \[ (x,y,z,d) = (-39.747837, -134.274144, -9413.624553, 0.0185173). \] The second solution exists because there are two solutions for the equation of 4 intersecting spheres.
As an alternative to using the multivariate Newton's method in the first step, the quadratic formula can be used to solve the system after careful algebraic manipulation of the equations listed in 4.38 of the reality check (also copied below). The steps can be found here and the corresponding code is available here . To begin, the second, third, and fourth equations are subtracted from the first equation to cancel all the quadratic terms. Then the equation is rewritten so the right side is 0. After doing this, the equations are reordered to the following form:
Ux x + Uy y + Uz z + Ud d + W = 0
where the vectors Ux, Uy, Uz, Ud, and W are the coefficients of the linear system. After writing this as a collection of linear equations, a careful determinant was used to simply x, y, and z to be represented as a function of d. Finally, with x, y, and z being functions of d, they are substituted into one of the equations below. The equation then reads F(d)=0 where F(d) is a messy quadratic polynomial. Lastly, the quadratic formula was applied to find the two possible roots of F(d).
Now set up a test of the conditioning of the GPS problem. Define satellite positions (Ai, Bi, Ci) from spherical coordinates (ρi, φi, θi) as
Ai = ρcos(φi)cos(θi)where ρ = 26750 km is fixed, while 0 ≤ φi ≤ π/2 and 0 ≤ θi ≤ 2π for i = 1,...,4 are chosen arbitrarily. The φ coordinate is restricted so that the four satellites are in the supper hemisphere. Set x = 0, y = 0, z = 6370, d = 0.0001, and calculate the corresponding satellite range and travel times.
Code can be found here.In our code, we use the following spherical coordinates: (26750, 0, 0), (26750, π/6, 2π/3), (26750, π/3, 4π/3), (26750, π/2, 2π).
We will define an error magnification factor specially tailored to the situation. The atomic clocks aboard the satellites are correct up to about 10 nanoseconds, or 10-8 seconds. Therefore, it is important to study the effect of changes in the transmission time of this magnitude. Let the backward, or input error be the input change in meters. At the speed of light, 10-8 seconds corresponds to roughly 3 meters. Let the forward, or output error be the change in position, caused by such a change in ti, also in meters. Then we can define the dimensionless
and the condition number of the problem to be the maximum error magnification factor for all small Δti.(say, 10-8 or less).
Change each ti defined in the foregoing by Δti = +10-8 or -10-8, not all the same. Denote the new solution of the equation by (x',y',z',d'), and compute the difference in position, on the basis of the error magnification factor.
With our given coordinates and guess solution, we found the following:
new solution(x',y',z',d') = (0.002832566559, 0.002026503684, 6370.004664877504, 0.000100015560)
variations of Δti: -1, -1, -1, 1
||Δx,Δy,Δz||∞ = 0.009329758378954
error magnification factor = 3.112072412093021
Now repeat Step 4 with a more tightly grouped set of satellites. Choose all φi within 5 percent of one another and all θi within 5 percent of one another. Solve with and without the same input error as in Step 4. Find the maximum position error and error magnification factor. Compare the conditioning of the GPS problem when the satellites are tightly or loosely bunched.
Code can be found here. All we had to do was to change two lines of code.In our code, we use the following spherical coordinates: (26750, 0, 0), (26750, π/40, 21π/40), (26750, π/20, 11π/20), (26750, 3π/40, 23π/40).
With our given coordinates and guess solution, we found the following:
new solution(x',y',z',d') = (-12.239979683903, 2.026503684, 6357.394422009528, 0.000102764332)
variations of Δti: 1, -1, 1, -1
||Δx,Δy,Δz||∞ = 16.846727318932608
error magnification factor = 5619.463355189746
As seen by the significantly larger error magnification factor, satellites bunched closely together are more likely to have errors in positioning than satellites that are more loosely distributed.
In parts 4 and 5 it was demonstrated that changes in time can increase the error of the approximate location calculated using the GPS equations. It is also possible for the error to be affected by the amount of satellites sending data to the reciever.
To determine the effect the amount of satellites has on the condition number and forward error of the GPS approximated location, we wrote programs to find the condition number and forward error for all amounts of satellites from 4 to 32. The first program written to find the data iterates through all satellite amounts 4 to 32. In the iteration, a program is used to find the condition number for each satellite amount by running a function that calculates the error magnification factor and forward error of a satellite amount 100 times. The function that finds the EMF and forward error does this by following the methods used in step 4 and 5. The computer programs used in this part are BestSatelliteNumber, SatelliteConditionNumber, SatelliteEMF, and the Newton's method equation from Part 1.
The same Newton's method used in part 1 is sufficient for the Gauss Newton method needed to find the coordinates from more than 4 satellites because the "backslash" command in Matlab uses the normal equations if the system has more functions than unknown values. In the Gauss Newton method, the normal equations iterated through to determine each new guess are \[ D^TDg=-D^TF \] where \(F\) is the vector of the function at the current guess, \(D\) is the Jacobian matrix of the function at the current guess, and \(g\) is the vector of unknown variables which will be the next guess.
Using the methods described in the two preceeding paragraphs, we found that using more than 4 satellites can reduce the condition number and forward error of the approximation. The results of the best satellite amounts are summarized in the following table:
| Satellite Amount | Characteristic | Value | 31 | Smallest Condition Number | 0.65372 | 32 | Lowest Median EMF | 0.29679 | 31 | Smallest Maximum Forward Error | 0.0019598 | 32 | Smallest Median Forward Error | 8.8975e-04 |