equation (1)
By: Austin Alderete, Theodore Faust, and Elisa Kim
Reference: Numerical Analysis by Timothy Sauer (p.238-241)
The global positioning system (GPS) consists of 24 satellites carrying atomic clocks, measuring time between the transmissions and the satellites in nanoseconds. The satellites are orbiting the earth at an altitude of 20,200km. Each satellite sends the signals, then the GPS receivers on earth obtain the information of the location from the transmission time and the position recieved from the signal. The distances between the receiver and the transmission can be calculated by multiplying it with speed of light, \(c\approx299792.458\) km/sec. However, because the time for the receiver is relatively inaccurate comparing to the atomic clock system, we need to add one more satellite to get the error in time, \(d\), in order to get more accurate position \((x,y,z\)). With the location of the satellite \((A_i,B_i,C_i\)), four equations to solve four variables \(x,y,z,d\) can be written as below:
The equations above can be rewritten as:
We can solve these equation using Multivariate Newtons Method, where we find the Jacobian Matrix, denoted Df, and solve the vector \(\vec v\) from \[ Df (\vec x_0)(\vec v) = f(\vec x_0) \] From the second equations, we got the Jacobian matrix as: \begin{matrix} 2(x-A_1) & 2(y-B_1) & 2(z-C_1) & 2c^2(t_1-d)\cr 2(x-A_2) & 2(y-B_2) & 2(z-C_2) & 2c^2(t_2-d)\cr 2(x-A_3) & 2(y-B_3) & 2(z-C_3) & 2c^2(t_3-d)\cr 2(x-A_4) & 2(y-B_4) & 2(z-C_4) & 2c^2(t_4-d) \end{matrix} and \(f(\vec x_0)\) can be written in matrix as: \begin{matrix} (x-A_1)^2 + (y-B_1)^2 + (z-C_1)^2 -[c(t_1-d)]^2\cr (x-A_2)^2 + (y-B_2)^2 + (z-C_2)^2 -[c(t_2-d)]^2\cr (x-A_3)^2 + (y-B_3)^2 + (z-C_3)^2 -[c(t_3-d)]^2\cr (x-A_4)^2 + (y-B_4)^2 + (z-C_4)^2 -[c(t_4-d)]^2 \end{matrix}
1. The four satellite positions are given:
(15600, 7540,20140), (18760, 2750, 18610), (17610, 14630, 13480), (19170, 610, 18390) in km.
The measured time intervals are 0.07074, 0.07220, 0.07690, 0.07242 in seconds.
The initial vector is \((x_0, y_0, z_0, d_0\)) = (0, 0, 6370, 0).
Note that we used n=10 iterations in the code. At the end of the iterations, we could get the desired result, \((x,y,z)=(-41.77271, -16.78919, 6370.0596)\) and \(d=-3.201566 \times 10^{-3}\).
We changed initial guess x0=[-600000; -600000; -600000; 0] and changed number of steps n=100
in the previous code.
And we got \((x,y,z)=(-39.747837348218, -134.274144360693, -9413.624553735819)\)
and \(d=0.185173047096\).
2. From the equation (2) above, we can expand equations and subtract last three equations from the first equation.
\begin{eqnarray}
x(-2A_1+2A_2) + y (-2B_1 +2B_2) +z(-2C_1-2C_2)-dc^2(-2t_1+2t_2)+({A_1}^2-{A_2}^2)+({B_1}^2-{B_2}^2)+ ({C_1}^2-{C_2}^2)-c^2({t_1}^2+{t_2}^2) &=& 0 \cr
x(-2A_1+2A_3) + y (-2B_1 +2B_3) +z(-2C_1-2C_3)-dc^2(-2t_1+2t_3)+({A_1}^2-{A_3}^2)+({B_1}^2-{B_3}^2)+ ({C_1}^2-{C_3}^2)-c^2({t_1}^2+{t_3}^2) &=& 0 \cr
x(-2A_1+2A_4) + y (-2B_1 +2B_4) +z(-2C_1-2C_4)-dc^2(-2t_1+2t_4)+({A_1}^2-{A_4}^2)+({B_1}^2-{B_4}^2)+ ({C_1}^2-{C_4}^2)-c^2({t_1}^2+{t_4}^2) &=& 0 \cr
\end{eqnarray}
Now, they are in the form of \(x\vec u_x + y\vec u_y +z\vec u_z + d\vec u_d+\vec w = 0\).
The formula for x in terms of d can be obtained from the equation given below:
\begin{eqnarray}
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] \cr
&=& \det [\vec u_y | \vec u_z | x \vec u_x ]+\det [\vec u_y | \vec u_z |y \vec u_y] +\det [\vec u_y | \vec u_z |z\vec u_z] +\det [\vec u_y | \vec u_z |d\vec u_d]+\det [\vec u_y | \vec u_z |\vec w] \cr
&=& x \det [\vec u_y | \vec u_z | \vec u_x ]+d \det[\vec u_y | \vec u_z |\vec u_d] +\det[\vec u_y | \vec u_z |\vec w] \cr
\end{eqnarray}
In this way, we can get all \(x,y,z\).
\begin{eqnarray}
x &=& (d \det([\vec u_y | \vec u_z |\vec u_d])+\det([\vec u_y | \vec u_z |\vec w]))/(-\det([\vec u_y | \vec u_z | \vec u_x]))\cr
y &=& (d \det([\vec u_z | \vec u_x |\vec u_d])+\det([\vec u_z | \vec u_x |\vec w]))/(-\det([\vec u_z | \vec u_x | \vec u_y]))\cr
z &=& (d \det([\vec u_x | \vec u_y |\vec u_d])+\det([\vec u_x | \vec u_y |\vec w]))/(-\det([\vec u_x | \vec u_y | \vec u_z]))\cr
\end{eqnarray}
The Matlab code can be found here.
NB! There is an alternative method of solving that uses Cramer's rule. The implementation of this method can be found within the code.
The far-from-earth solution to the equations: \((x,y,z)=(-39.747837348143 -134.274144360667 -9413.624553735755)\) in km.
\(d=0.000185173047096\) seconds.
4. In this probelm, we defined satellite positions \((A_i,B_i,C_i\)) from spherical coordinates \((\rho,\phi_i,\theta_i)\),
as
\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 we set \(\rho=42256\). We restricted \(\phi_i\) to be between 0 and \(\frac{\pi}{2}\).
We set \(\phi_i\) to be \(\phi_1=0, \phi_2=\frac{\pi}{8}, \phi_3=\frac{\pi}{4}, \phi_4=\frac{3\pi}{8} \)
and \(\theta_i\) to be \(\theta_1=0, \theta_2=\frac{\pi}{2}, \theta_3=\pi, \theta_4=\frac{3\pi}{2}\).
After setting \((x, y, z, d) = (0, 0, 6370, 0.0001)\), we calculated the satellite ranges and travle times with equations followed:
\begin{eqnarray}
R_i &=&\sqrt{{A_i}^2+{B_i}^2+{(C_i-6370)}^2}\\
t_i &=& d+\frac {R_i}c \\
\text{error magnification factor} &=& \frac {\|(\Delta x, \Delta y, \Delta z)\|_\infty}{c\|\Delta t_1,..., \Delta t_m \|_\infty} \\
\end{eqnarray}
The Results of EMFs(error magnification factors) are found by running the code:
\begin{array}{|c|l|}
\hline
(t_1,t_2,t_3,t_4)^* & EMF \\ \hline
(+,+,+,+) & 9.101243318531999\times 10^{-9} \\ \hline
(+,+,+,-) & 3.324033147145058 \\ \hline
(+,+,-,+) & 2.699141681387724 \\ \hline
(+,-,+,+) & 2.116970337014859 \\ \hline
(-,+,+,+) & 2.145221783231499 \\ \hline
(+,+,-,-) & 1.126352492950904 \\ \hline
(+,-,+,-) & 4.844362543876775 \\ \hline
(+,-,-,+) & 1.636669005810785 \\ \hline
(-,+,+,-) & 1.636668881423539 \\ \hline
(-,+,-,+) & 4.844364012514071 \\ \hline
(-,-,+,+) & 1.126352499189273 \\ \hline
(+,-,-,-) & 2.145221668555834 \\ \hline
(-,+,-,-) & 2.116970795584959 \\ \hline
(-,-,+,-) & 2.699141422305664 \\ \hline
(-,-,-,+) & 3.324034073651628 \\ \hline
(-,-,-,-) & 1.061811720495400\times 10^{-8} \\ \hline
\end{array}
* Postive sign means \(\Delta t_i = +10^{-8}\) and negative sign means \(\Delta t_i = -10^{-8}\), where \(i\)=1,2,3,4.
The condition number (the maximum error magnification factor) is approximately 4.844.
The maximum position error is found in (-,+,-,+) case, which is at (0.006279970976591, 0.006706432273050, 0.014523037947583) in km.
Therefore, (6.279970976591, 6.706432273050, 14.523037947583) in meter.
5. In order to make the satellites close to each other, we set \(\phi_i\) to be \(\phi_1=\frac{\pi}{4}, \phi_2=\frac{1.01\pi}{4}, \phi_3=\frac{1.02\pi}{4}, \phi_4=\frac{0.99\pi}{4} \) and \(\theta_i\) to be \(\theta_1=\pi, \theta_2=1.01\pi, \theta_3=1.02\pi, \theta_4=0.99\pi \).
The Results of EMFs are found by running the code: \begin{array}{|c|l|} \hline (t_1,t_2,t_3,t_4) & EMF \\ \hline (+,+,+,+) & 0.001293172910022 \\ \hline (+,+,+,-) & 2.436716303587290\times 10^5 \\ \hline (+,+,-,+) & 2.404587674552902\times 10^5 \\ \hline (+,-,+,+) & 7.546248790109969\times 10^5 \\ \hline (-,+,+,+) & 7.002846204809517\times 10^5 \\ \hline (+,+,-,-) & 2.528868894377193\times 10^3 \\ \hline (+,-,+,-) & 1.017901473347243\times 10^6 \\ \hline (+,-,-,+) & 4.956863410637114\times 10^5 \\ \hline (-,+,+,-) & 4.727176717485383\times 10^5 \\ \hline (-,+,-,+) & 9.255674489655442\times 10^5 \\ \hline (-,-,+,+) & 2.528712842781894\times 10^3 \\ \hline (+,-,-,-) & 7.519785962566979\times 10^5 \\ \hline (-,+,-,-) & 7.027282170913572\times 10^5 \\ \hline (-,-,+,-) & 2.462357419081989\times 10^5 \\ \hline (-,-,-,+) & 2.379617217600075\times 10^5 \\ \hline (-,-,-,-) & 3.605821590369192\times 10^{-5} \\ \hline \end{array}
The condition number is approximately 1.0179\(\times 10^6\).
The maximum position error achieved in the (+,-,+,-) case, which was (195.475089709398, 802.719725023311, 3051.591846965914) in km.
Therefore, (195475.089709398, 802719.725023311, 3051591.846965914) in meter.
When the GPS satellites are close together, the condition number was more than \(10^5\) times larger than
when they were spread apart. This led to errors of 3000 km rather than 0.014 km.
6. In this problem, We tested whether the GPS error and condition number can be reduced by adding more satellites. We ran the code for different number of satellites from 1 to 10 to see which case gives us the best result. since we have 4 variables in the equations, we can't have less than four satellites in order to solve equations to get the location. \begin{array}{|c|c|} \hline \text{Number of Satellites}\ & \text{Maximum EMF} \ \\ \hline 1 & \text{N/A} \ \\ \hline 2 & \text{N/A} \ \\ \hline 3 & \text{N/A} \ \\ \hline 4 & 4.844364021918690 \\ \hline 5 & 4.292302631019335 \\ \hline 6 & 3.674773851463355 \\ \hline 7 & 3.883475386383558 \\ \hline 8 & 3.808996695018498 \\ \hline 9 & 3.532017744061210 \\ \hline 10 & 3.597401154938986 \\ \hline \end{array} As our result, we have the lowest error when we had 9 setellites with the maximum EMF= 3.532.