Project 2: GPS, Conditioning, and Nonlinear Least Squares

Team members: Brendan Gramp, Michael Belete, Faysal Shaikh

Contents

Background

The global positioning system (GPS) consists of 24 satellites carrying atomic clocks, orbiting the earth at an altitude of 20,200 km. Four satellites in each of six planes, slanted at $55^\circ$ degrees with respect to the poles, make two revolutions per day. At any time, from any point on earth, five to eight satellites are in the direct line of sight. Each satellite has a simple mission: to transmit carefully synchronized signals from predetermined positions in space, to be picked up by GPS receivers on earth. The receivers use the information, with some mathematics, to determine accurate $(x,y,z)$ coordinates of the receiver.

View Reality Check 4.

Part 1-- Obtain numerical solutions for the system.

Solve the system by using Multivariate Newtons Method. Find both receiver positions $(x,y,z)$ near to and far from earth, with time correction $d$ values, for known, simultaneous satellites [in the form of $(x,y,z,d)$]: $(15600,7540,20140,0.07074)$, $(18760,2750,18610,0.07220)$, $(17610,14630,13480,0.07690)$, $(19170,610,18390,0.07242)$. Set the initial vector to be $(0,0,6370,0)$.

A function, gps.m, was created to perform multivariate Newtons method given an input of 4 column vectors of the form $(x,y,z,d)$ and $n$, the number of iterations to follow using multivariate Newtons method. A 5th column vector containing an initial guess, $(x_0,y_0,z_0,d_0)$, for multivariate Newtons method is hardcoded-in.

View the full function, gps.m.

Using our initial guess vector described above (and hardcoded-in), we first initialize the vectors of the four satellites given above. We then initialize the counter for multivariate Newtons method (e.g., we wanted 10 iterations of multivariate Newtons method here, so we set n=10). Finally, we run the function with its appropriate inputs and examine the answer.

%format long to view more decimal places
format long

%initialize vectors
v1=[15600; 7540; 20140; 0.07074];
v2=[18760; 2750; 18610; 0.07220];
v3=[17610; 14630; 13480; 0.07690];
v4=[19170; 610; 18390; 0.07242];

%initialize multivariate Newtons method 
counter
n=10;

%run multivariate Newtons using v1,v2,v3,v4,n
%and initial guess (0,0,6370,0)[hardcoded-in]
near_solution=gps(v1,v2,v3,v4,n)

near_solution =

   1.0e+03 *

  -0.041772709570857
  -0.016789194106537
   6.370059559223330
  -0.000003201565830

While the above is promising (a solution is indeed produced), we had a difficult time attempting to produce the supposed second solution of the system. As the initial guess was harcoded-in for the above gps.m function, it would be tedious to show each newly-defined function with a different initial guess.

For description's sake, several different initial guesses for multivariate Newtons method were provided, yet the system either diverged to infinity, or converged to the solution shown above.

Part 2-- Obtain exact solutions for the system.

Write a MATLAB program to solve the problem in step 1 via the quadratic formula. This was completed following the steps described in the revised version of Step 2.

View revised Step 2.

A function, gpsquadratic.m, was defined to return both exact solutions of the system when converted into a quadratic equation. The function requires input of 4 column vectors of the form $(x,y,z,d)$ [just as gps.m], but since we are finding exact solutions, we no longer require the iterator $n$ [unlike for gps.m].

View the full function, gpsquadratic.m.

Here we input the previously-utilized input vectors to obtain both solutions of the system:

%solve quadratic 
for exact solutions of system
[first_sol,second_sol]=gpsquadratic(v1,v2,v3,v4)

first_sol =

   1.0e+03 *

  -0.041772709570836
  -0.016789194106526
   6.370059559223344
  -0.000003201565830


second_sol =

   1.0e+03 *

  -0.039747837348165
  -0.134274144360667
  -9.413624553735767
   0.000185173047096

We seemingly found the value of first_sol in Part 1, but had issues with having multivariate Newtons method converge to second_sol in Part 1. This was because our Newtons method iterator used equations 4.37 (the square-root forms) and for Part 2 our quadratic is from equations 4.38 (the squared forms).

Part 3-- Use MATLAB's Symbolic Toolbox.

This part was omitted for this project.

Part 4-- Perform a test of condtioning of the GPS problem.

Define satellite positions $(A_i,B_i,C_i)$ from spherical coordinates $(\rho,\phi_i,\theta_i)$ as $A_i=\rho\cos\phi_i\cos\theta_i$, $B_i=\rho\cos\phi_i\sin\theta_i$, and $C_i=\rho\sin\phi_i$, where $\rho=26570$ kilometers is fixed, while $0\leq\phi_i\leq\frac{\pi}{2}$ and $0\leq\theta_i\leq2\pi$ for $i=1,...,4$ are chosen arbitrarily. The $\phi$ coordinate is restricted so that the four satellites are in the upper hemisphere. Set $x=0$, $y=0$, $z=6370$, $d=0.0001$, and calculate the corresponding satellite ranges $R_i=\sqrt{A^2_i+B^2_i+(C_i-6370)^2}$ and travel times $t_i=d+\frac{R_i}{c}$.

We will define an error magnification factor (emf) specially tailored to the situation. The atomic clocks aboard the satellites are correct up to about 10 nanoseconds, or $10^{-8}$ second. 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, $\Delta t_i=10^{-8}$ second corresponds to $10^{-8}c\approx3$ meters. Let the forward, or output error be the change in position $||(\Delta x,\Delta y,\Delta z)||_\infty$, caused by such a change in $t_i$, also in meters. Then we can define the dimensionless $emf=\frac{||(\Delta x,\Delta y,\Delta z)||_\infty}{c||(\Delta t_1,..., \Delta t_m)||_\infty}$, and the condition number of the problem to be the maximum emf for all small $\Delta t_i$ (say, $10^{-8}$ or less).

Change each $t_i$ defined in the foregoing by $\Delta t_i=+10^{-8}$ or $-10^{-8}$, not all the same. Denote the new solution of the equations by $(x,y,z,d)$ and compute the difference in position $||(\Delta x,\Delta y,\Delta z)||_\infty$ and the emf. Try different variations of the $\Delta t_i$'s. What is the maximum position error found, in meters? Estimate the condition number of the problem, on the basis of the emf's you have computed.

We first began by defining a function, gpsspherical.m, for use in spherical coordinates. Similar to previous gps.m and gpsquadratic.m functions, gpsspherical.m requires an input of 4 vectors (this time in the form of $(\phi_i,\theta_i)$, since the problem states that $\rho=26570$ kilometers is fixed for all 4 satellites (all satellites exhibit equal radii: their distances from the center of the earth). The function essentially serves to convert spherical coordinates given into a set of rectangular coordinates, with the addition of error given by $\pm\Delta t_i$. The addition or subtraction of $\Delta t_i$ was decided and hardcoded-in for each attempt. Thus, the above-shown code only describes one distribution of $+\Delta t_i$ and $-\Delta t_i$. The final results of the function are to compute both the difference in position, $||(\Delta x,\Delta y,\Delta z)||_\infty$, and the emf, and $\frac{||(\Delta x,\Delta y,\Delta z)||_\infty}{c\cdot10^{-8}}$.

View the full function, gpsspherical.m

We begin by initializing the values in spherical coordinates describing the locations of the satellites.

%initialization:
s1 = [pi/6; pi/6];
s2 = [pi/3; pi/12];
s3 = [3*pi/2; pi/4];
s4 = [5*pi/6; pi/9];

We also visualize our satellites to show they are in an unclustered configuration around the earth. We do this by first creating a function, gpsplot.m, which allows us to visualize 4 satellites around the earth.

View the full function, gpsplot.m

%visualization of satellites:
gpsplot(s1,s2,s3,s4)

Then, we run the function, in order to determine values for the error in position and the emf.

%run:
unclust=gpsspherical(s1,s2,s3,s4)

unclust =

  19.690451877553240
   0.059030489675024

Here, we notice that our positional error is approx 0.059km, which converts to 59m. This is a modest error, but likely wouldn't result in many practical issues for GPS purposes.

Part 5-- Repeat Step 4 with a more tightly grouped set of satellites.

Choose all $\phi_i$ within 5% 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.

We utilize the same gpsspherical.m function described above, but intialize the function with different values for the spherical coordinates.

%initialize, with all coordinates 
within 5% of one another:
s1 = [pi/6; pi/6];
s2 = [(pi/6)*1.05; (pi/6)*.95];
s3 = [(pi/6)*1.025; (pi/6)*.975];
s4 = [(pi/6)*1.0375; (pi/6)*1.03];

As in Part 4, we would like to visualize the configuration of the satellites around the earth. We do so with the same gpsplot.m function used earlier.

%visualization of satellites:
gpsplot(s1,s2,s3,s4)

And we then re-run the same function to obtain our new answers for clustered satellites.

%run:
clust=gpsspherical(s1,s2,s3,s4)

clust =

   1.0e+03 *

   9.280243248441554
   0.027821469342882

Positional error is now approx. 27.8km, which would likely be problematic for the vast majority of GPS uses.

Part 6-- Decide if adding satellites helps conditioning/error.

Return to the unbunched satellite configuration of Step 4, and add more satellites; any number besides 4 additional satellites is acceptable. (At all times and at every position on earth, 5 to 12 GPS satellites are visible.) Design a Gauss-Newton iteration to solve the least squares system of eight equations in four variables $(x,y,z,d)$. What is a good initial vector? Find the maximum GPS position error, and estimate the condition number. Summarize your results from four unbunched, four bunched, and Part 6's unbunched satellites. What configuration is best, and what is the maximum GPS error, in meters, that you should expect solely on the basis of satellite signals?

We decided to investigate adding 5 satellites to the current system, making the total number of satellites equal to 9. We redefine satellites we have used from Part 4, as well as the positions we have defined for our new satellites.

%satellite 
position assignment
s1 = [pi/6; pi/6];
s2 = [pi/3; pi/12];
s3 = [3*pi/2; pi/4];
s4 = [5*pi/6; pi/9];
s5 = [1.05*pi; 41*pi/96];
s6 = [pi/5; pi/3.5];
s7 = [2.121*pi/2; 3*pi/8];
s8 = [5.5*pi/6; pi/3.3];
s9 = [pi/5; pi/2.1];

In order to visualize the problem with adding 5 satellites (making the total number of satellites equal to 9), we define a new function, gpsplot2.m.

View the full function, gpsplot2.m

We then attempt to visualize the configuration of all 9 satellites around the earth.

%visualization of 
satellites:
gpsplot2(s1,s2,s3,s4,s5,s6,s7,s8,s9)

We redefine a new Gauss-Newton method in order for us to iterate towards a solution. We do so by defining a new function, gps3.m. It should be noted that our chosen initial vector, hardcoded-in, is (0,0,6370,0). This seemed to work reasonably-well for our purposes.

View the full function, gps3.m

However, our current intputs are in the form of spherical coordinates. [|gps3.m| requires rectangular coordinates.] We define a new function, gpsspherical2.m, for the spherical coordinates conditioning problem answered in Part 4 and Part 5. [This function, in a much similar manner to gpsspherical.m, also converts spherical coordinates into rectangular coordinates, and then runs them through the iterative numerical method, which is gps3.m in this case.]

View the full function, gpssherical2.m

And we attempt to run the same conditioning problem, using the spherical coordinates of our 4 original and 5 newly added satellites.

%run with 9 satellites
part6=gpsspherical2(s1,s2,s3,s4,s5,s6,s7,s8,s9)

part6 =

   3.292673355575251
   0.009871186386590

This run shows us that with these 9 total satellites, we would have a positional error of approximately 0.010km, which converts to 10m. Our test of 9, unbunched satellites seems to result in a better conditioning of the system, and thus a smaller emf and positional error.


Published with MATLAB® R2017a