Project 2
GPS, Conditioning, and Nonlinear Least Squares
Lucas Bouck, Patrick Bishop, and Paul McNulty
The problem here was create a method to obtain an accurate global position, and explore
what aspects could affect the accuracy of those positions. GPS satellites send signals
with information such as location and time. In order for a (semi) unique position to be
obtained, the signals from four such satellites must be received by a positioning device,
three for the spatial (x,y,z) unknowns and one to make sure the temporal accuracy is
decent. These signals can be used to construct a system of nonlinear equations which can
be solved to give a location of the device using Multivariate Newton's Method or Nonlinear
Least Squares.
Part 1
In the first part we were asked to find the location of the device given the location of 4
satellites using multivariate Newton's method. The relevant equations used are:
\begin{eqnarray*}
(x-A_1)^2+(y-B_1)^2+(z-C_1)^2 &=& c^2(t_1-d)^2 \\
(x-A_2)^2+(y-B_2)^2+(z-C_2)^2 &=& c^2(t_2-d)^2 \\
(x-A_3)^2+(y-B_3)^2+(z-C_3)^2 &=& c^2(t_3-d)^2 \\
(x-A_4)^2+(y-B_4)^2+(z-C_4)^2 &=& c^2(t_4-d)^2
\end{eqnarray*}
Where each \( (A_i,B_i,C_i)\) are the spatial coordinates for each satellite and each
\(t_i\) is the respective temporal component. Our goal hear is to solve for \((x,y,z,d)\),
i.e. the location of the device.
Multivariate Newton requires an initial guess, which we pass in as a 4 dimensional vector
\(x_0\). Each vector \(x\) is passed into the function F (the four equations above), and
we must also compute the Jacobian of F, which is the matrix J, and is given by:
$$J=\left[\begin{matrix}\frac{\partial f_1}{\partial x} & \frac{\partial f_1}{\partial y} & \frac{\partial
f_1}{\partial
z} & \frac{\partial f_1}{\partial d}\\
\frac{\partial f_2}{\partial x} & \frac{\partial f_2}{\partial y} & \frac{\partial f_2}{\partial
z} & \frac{\partial f_2}{\partial d}\\
\frac{\partial f_3}{\partial x} & \frac{\partial f_3}{\partial y} & \frac{\partial f_3}{\partial
z} & \frac{\partial f_3}{\partial d}\\
\frac{\partial f_4}{\partial x} & \frac{\partial f_4}{\partial y} & \frac{\partial f_4}{\partial
z} & \frac{\partial f_4}{\partial d}\\
\end{matrix}\right]$$
where each component in J is a partial derivative of F with respect to a variable. From
there, it is a simple algorithm:
\begin{eqnarray*}
J(x_n)v &=& F(x_n)\\
x_{n+1} &=& x_{n}-v
\end{eqnarray*}
The matrix equation is then solved using gaussian elimination. The code can be seen here:
newtonmgps1.m.
Using an initial guess of (0, 0, 6370, 0) and satellite coordinates (15600, 7540, 20140, 0.07074),
(18760, 2750, 18610, 0.07220), (17610, 14630, 13480, 0.07690), and (19170, 610, 18390, 0.07242)
with the spatial coordinates measured in kilometers and the temporal component measured in
seconds, our code found the solution, \((x,y,z,d)=(-41.772709570875,-16.789194106543,6370.059559223320,
-3.201565830*10^{-3})\), which is what we wanted.
Part 2
Part 2 asked us to use a different method to solve for the same thing as Part 1. Using
determinants and the same function F from part one, we were able to algebraically isolate
one variable at a time which we could solve for. By subtracting one of the bottom three
equations from the first, we were able to come up with three new equations, resulting in a
vector equations $\vec{u_x}x+\vec{u_y}y+\vec{u_z}z+\vec{w}=0$. We could then create an
equation for x (or y,z) in terms of d by allowing $0=\det[\vec{u_y}|\vec{u_z}|\vec{u_x}x+\vec{u_y}y+\vec{u_z}z+\vec{w}]$.
Once one variable is found, we could use the equation $$\sqrt{(x-A_1)^2+(y-B_1)^2+(z-C_1)^2}-c(t_1-d)=0$$
and the bisection method to solve for d, and thus for x, y, and z as well.
Using the same test case as in Part 1, we were able to use the above method to find the
correct solution. The relevant code can be found here: part2.m
and f2.m.
Part 4
Part 4 asked us to look at the conditioning of the GPS problem. Each satellite can be
represented in spherical coordinates given by:
\begin{eqnarray*}
A_i &=& \rho \cos(\phi_i) \cos(\theta_i)\\
B_i &=& \rho \cos(\phi_i) \sin(\theta_i)\\
C_i &=& \rho\sin(\phi_i)
\end{eqnarray*}
Where now $\rho$ is the radius from the center of the earth to the satellite, $\phi$ gives
the altitude as measured from the north pole, and $\theta$ is the azimuthal angle. For
this problem we will restrict
$\rho=26570$, $0\leq \phi_i\leq \frac{\pi}{2}$, and $0\leq \theta_i\leq 2\pi$ for all $i=1,2,3,4$.
To test the conditioning, we picked the known location from Part 1, and then picked
satellites that were fairly evenly spaced out but still within our restrictions. We then
had to calculate the range of each satellite to our location and the time it takes for the
signal to reach our locations. The relevant equations for this are:
\begin{eqnarray*}
R_i &=& \sqrt{A_i^2+B_i^2+(C_i-6370)^2}\\
t_i &=& d+R_i/c
\end{eqnarray*}
The mostly likely error will be found in the temporal component, so that is where we will
base our conditioning. We will change all of the satellite times by a combination of
$\pm1*10^{-8}$, which $10^{-8}$ is approximately the precision the clocks on the
satellites are capable of. Once we have the calculated postion, we will compare it to the
original positions using the error vector $(\Delta x,\Delta y,\Delta z)$. The error will
be measured using the infinit norm, such that max error $=||(\Delta x,\Delta y,\Delta z)||_{\infty}$.
This will be reported in meters. We will then find the magnification error given by
$$\frac{||(\Delta x,\Delta y,\Delta z)||_{\infty}}{c||(\Delta t_1, \Delta t_2, \Delta t_3, \Delta t_4)||_{\infty}}$$
The error magnification will give us a ratio of the max position error to the distance that
light travels in the amount of time that we are changing it by. We ran eight runs with
different combinations of time errors. Code for this part can be found here:
part4.m. Our results are given below:
| Trial |
1 |
2 |
3 |
4 |
5 |
6 |
7 |
8 |
Maximum |
| Max Position Error (m) |
12.1807 |
25.2583 |
25.4063 |
38.4840 |
12.1807 |
25.2583 |
25.4063 |
38.4839 |
38.4839 |
EMF |
4.063 |
8.4253 |
8.4746 |
12.8369 |
4.0630 |
8.4253 |
8.4746 |
12.8369 |
12.8369 |
From which we can conclude that the condition number for this problem is about 12.8369.
Part 5
For Part 5, we again want to find the error and conditioning number, but this time our
satellites are bunched up in the sky, and so the problem is ill-conditioned. All satellites
were within five perent of each other, and the code can be found here:
part5.m. The results are given below:
| Trial |
1 |
2 |
3 |
4 |
5 |
6 |
7 |
8 |
Maximum |
| Max Position Error (m) |
5510 |
8262 |
14,533 |
743 |
5514 |
8270 |
14,508 |
743 |
14,533 |
EMF |
1838 |
2756 |
4848 |
248 |
1839 |
2759 |
4839 |
248 |
4848 |
From which we can conclude a condition number of 4848.
Part 6
For part 6, we were to use eight satellites to get out position. The method was the same
as in part 4, and this time we did sixteen trials in order to account for more combinations
of satellite errors. However, since we have many more unknowns, Newton's Method is not
sufficient for solving this problem. Instead, we used Nonlinear Least Squares. The algorithm
is analogous to Newton's method, except now the matrix equation is $$J^{T}Jv=J^{T}F$$.
Where $J^T$ is now the transpose of the Jacobian matrix. The code for part six can be found
here: newtonmgps6.m, and the positions are given by:
part6.m. The results are given below:
| Trial |
1 |
2 |
3 |
4 |
5 |
6 |
7 |
8 |
9 |
10 |
11 |
12 |
13 |
14 |
15 |
16 |
Maximum |
| Max Position Error (m) |
2.7839 |
2.2958 |
2.8388 |
2.3925 |
0.6332 |
0.6332 |
0.6332 |
0.6332 |
2.7839 |
2.2958 |
2.8388 |
2.3925 |
0.6332 |
0.6332 |
0.6332 |
0.6332 |
2.8388 |
EMF |
0.9286 |
0.7658 |
0.9470 |
0.7981 |
0.2112 |
0.2112 |
0.2112 |
0.2112 |
0.9286 |
0.7658 |
0.9470 |
0.7981 |
0.2112 |
0.2112 |
0.2112 |
0.2112 |
0.9470 |
Thus the condition number for part 6 was .9470, and therefore much better. The maximum
position error was only 2.8388, which is very accurate. Thus having 8 satellites is better
than only having four.