Computer Algebra Solution of GPS N-Points Problem

Béla  Palánc(z^1), Joseph L. Awang(e^2) and Erik Grafaren(d^3)

Abstract

In this contribution we illustrate different methods available for solving Global Positionig System N-points problem  using built-in functions of  Computer Algebra Systems (CAS). Direct and indirect minimization methods, Gauss-Jacobi combinatorial algorithm  with lexicographic and reduced Groebner basis, as well as with Dixon resultant and Sturmfels` formulation are demonstrated. Numeric and symbolic solutions are given for realistic situations. Using Mathematica 5.2 System, numeric and symbolic solution needs less computation time than that of iterative direct numerical method.

Keywords: GPS, CAS, Positioning, Algebra

1.  Introduction

In the present world, the art of positioning is changing rapidly with more Global Navigation Satellite System (GNSS), e.g., GALILEO (The European equivalent of GPS) being developed.  GNSS positioning is a modern surveying tool designed to give three dimensional coordinates in space. It has become a house hold name and its use traverse from military to weather parameter monitoring. Whereas access to GPS receivers used to be a preserve of those who could afford, wrist watches, mobile phones and other portable devices are fast being fitted with GPS receivers. A push of the button is sufficient to locate the user’s position. Behind the push of the button underlines complex mathematics which has to be manipulated to offer solution.
For users who would wish to go a step further than pushing the button to solve for the receiver position from the raw GNSS data, for purposes such as ambiguity resolution or nonlinear error modelling, we have here provide Mathematica 5.2 source codes (notebook) which could assist them in solving GNSS related problems. The codes solve numerical, semi-symbolic and symbolic multidimensional N-point GNSS problems. Though GPS is used as an example, the codes are however applicable to all GNSS related problems. It is hoped that the materials herein will provide you, the user, with information that will help you understand the behind the scene manipulation of GPS data.

2.  Ranging by GPS

Throughout history, position determination has been one of the most important tasks of mountaineers, pilots, sailor, civil engineers etc. In modern times, Global Positioning System (GPS) employing Global Navigation Satellite Systems (GNSS) provide an ultimate method to accomplish this task. If one has a hand held GPS receiver, the receiver measures the travel time of the signal transmitted from the satellites. Then this distance can be computed by multiplying the measured time by the speed of light in vacuum. The distance of the receiver from the i-th satellite, (Subscript[d, i]) is related to the unknown position of the receiver, {X,Y, Z}

\(d\_i = \@\(\((a\_i - X)\)\^2 + \((b\_i - Y)\)\^2 + \((c\_i - Z)\)\^2\) + \ \ ξ\)

where {(Subscript[a, i]),(Subscript[b, i]),(Subscript[c, i])},i =0.  .. 3 are the coordinates of the i-th satellite.

The distance is influenced also by the satellite and receiver' clock biases. The satellite clock biases can be modelled while the receiver' clock biases have to be considered as an unknown variable, ξ. This means, we have four unknowns, consequently we need four satellite signals as minimum observation. Let us employ (Subscript[x, 1]), (Subscript[x, 2]), (Subscript[x, 3]) and (Subscript[x, 4]) variables for the four unknowns, {X,Y, Z, ξ,} then our equations are, (Subscript[e, i]) = 0, i = 1...4, where

e1 = (((x1 - (Subscript[a, 0])))^2) + (((x2 - (Subscript[b, 0])))^2) + (((x3 - (Subscript[c, 0])))^2) - (((x4 - (Subscript[d, 0])))^2);

e2 = (((x1 - (Subscript[a, 1])))^2) + (((x2 - (Subscript[b, 1])))^2) + (((x3 - (Subscript[c, 1])))^2) - (((x4 - (Subscript[d, 1])))^2);

e3 = (((x1 - (Subscript[a, 2])))^2) + (((x2 - (Subscript[b, 2])))^2) + (((x3 - (Subscript[c, 2])))^2) - (((x4 - (Subscript[d, 2])))^2);

e4 = (((x1 - (Subscript[a, 3])))^2) + (((x2 - (Subscript[b, 3])))^2) + (((x3 - (Subscript[c, 3])))^2) - (((x4 - (Subscript[d, 3])))^2);

Let us suppose that the observation data are, see Grafarend and Shan [14],

data = {(Subscript[a, 0]) -> 1.483230866 (10^7), 
          (Subscript[a, 1]) -> -1.579985405 (10^7), 
          (Subscript[a, 2]) -> 1.98481891 (10^6), 
          (Subscript[a, 3]) -> -1.248027319 (10^7), 
(Subscript[b, 0]) -> -2.046671589 (10^7), (Subscript[b, 1]) -> -1.330112917 (10^7), 
          (Subscript[b, 2]) -> -1.186767296 (10^7), 
(Subscript[b, 3]) -> -2.338256053 (10^7), 
(Subscript[c, 0]) -> -7.42863475 (10^6), (Subscript[c, 1]) -> 1.713383824 (10^7), 
          (Subscript[c, 2]) -> 2.371692013 (10^7), 
(Subscript[c, 3]) -> 
            3.27847268 (10^6),  
(Subscript[d, 0]) -> 
            2.4310764064 (10^7), (Subscript[d, 1]) -> 2.2914600784 (10^7), 
          (Subscript[d, 2]) -> 2.0628809405 (10^7), 
(Subscript[d, 3]) -> 
            2.3422377972 (10^7)};

First, we introduce  different solution methods for this 4-points problem.

3.  Solutions of  the 4-points problem

3.1  Direct numerical solution

The most simple, direct solution is the direct numerical solution. Employing the built-in function of Mathematica, NSolve, we get

sol = 
      NSolve[Map[(# /. data) == 0 &, {e1, e2, e3, e4}], {x1, x2, 
            x3, x4}] // Timing

The results indicate that the solution is not unique. In order to decide on the admissible solution from the pair above, we compute the radial distance from the center of the Earth, the norm of the two positional vectors (Subscript[{(Subscript[x, 1]),(Subscript[x, 2]), (Subscript[x, 3])}, j]), j = 1, 2, namely,

Map[
      Norm[Drop[#, -1]] &, {x1, x2, x3, x4} /. sol[[2]]]

If the receiver coordinates are in the global reference frame, the norm of the positional vector of the receiver station will approximate the value of the Earth's radius. That is why the second solution is admissible,

SetPrecision[sol[[2, 2]], 10]

In order to check the result, let us  backsubstitute these values into the original equations,

{e1, e2, e3, e4} /. data /. sol[[2, 2]]

which means, that the absolute value of the highest residual is less than 20 centimeter. However, in this case, we desire small error rather than small residuals. The reader can easily verify, that computing the solution with the significance-precision of the input data, with 10 digits, instead of the default machine precision 16 digits, these residuals will be zero!

In order to check the error of the solution, one may repeat this computation with a higher precision,  i.e. with 30 digits.

data30 = 
        Map[#[[1]] -> SetPrecision[#[[2]], 30] &, data];

sol30 = 
        NSolve[Map[(# /. data30) == 0 &, {e1, e2, e3, e4}], {x1, 
            x2, x3, x4}, 30];

sol30[[2]]

SetPrecision[sol30[[2]], 10]

There is  no difference between the two solutions, however now the residuals disappeared,

{e1, e2, e3, e4} /. data30 /. sol30[[2]]

3.2  Symbolic solution

First, this equation  system of polynomials, will be transformed into  a system of linear equations and a quadratic equation. Let us expand and multiply by minus one, and arranged the original equations,

eqsL = {e1L, e2L, e3L, e4L} = 
        Map[-Sort[Expand[#]] &, {e1, e2, e3, e4}]

Subtruct the fourth equation from the other three ones,

q = 
      Table[eqsL[[i]] - eqsL[[4]], {i, 1, 3}] // Simplify

This is a system of three linear equations, and they can be written as,

(Subscript[g, 1]) = (Subscript[a, 0, 3]) x1 + (Subscript[b, 0, 3]) x2 + (Subscript[c, 0, 3]) 
            x3 + (Subscript[d, 3, 0]) x4 + (Subscript[e, 0, 3]);

(Subscript[g, 2]) = (Subscript[g, 1]) /. {1 -> 2, 0 -> 1}

(Subscript[g, 3]) = (Subscript[g, 1]) /. {1 -> 3, 0 -> 2}

The coefficients {(Subscript[a, i, 3]), (Subscript[b, i, 3]), 
          (Subscript[c, i, 3]), (Subscript[d, 3, i]), (Subscript[e, i, 3])}, i = 0.  ..  2, can be determined  as,

coeffs0 = 
        Table[{Coefficient[q[[i + 1]], x1], 
              Coefficient[q[[i + 1]], x2], 
              Coefficient[q[[i + 1]], x3], 
              Coefficient[q[[i + 1]], x4]} // Factor, {i, 0, 
            2}];

which are the coefficients of the variables {(Subscript[x, 1]),(Subscript[x, 2]), (Subscript[x, 3]), (Subscript[x, 4])}. The constant part is,

coeffs1 = 
        Table[{q[[i]] - coeffs0[[i]] . {x1, x2, x3, x4} // 
              Simplify}, {i, 1, 3}];

Therefore, all of the coefficients are,

coeffs = 
      Table[Union[coeffs0[[i]], coeffs1[[i]]], {i, 1, 3}]

Let us assign these coefficients to the linear system,

coeffsn = 
      Flatten[Table[
          Inner[#1 -> #2 &, {(Subscript[a, i, 3]), (Subscript[b, i, 3]), (Subscript[c, i, 3]), 
              (Subscript[d, 3, i]), (Subscript[e, i, 3])}, coeffs[[i + 1]], List], {i, 0, 
            2}]]

In addition, we take one of the nonlinear equations, let say the fourth one,

(Subscript[e, 4]) = (((x1 - (Subscript[a, 3])))^2) + (((x2 - (Subscript[b, 3])))^2) + (((x3 - (Subscript[c, 3])))^2) - (((x4 - (Subscript[d, 3])))^2);

Now, we shall solve the linear system for the variables {(Subscript[x, 1]),(Subscript[x, 2]), (Subscript[x, 3])}, with (Subscript[x, 4]) as parameter. It means, the relations (Subscript[x, 1])= g((Subscript[x, 4])), (Subscript[x, 2])= g((Subscript[x, 4])) and (Subscript[x, 3])= g((Subscript[x, 4])) will be computed. To do that, different methods can be employed.

3.2.1 Sturmfels' method

The Sturmfels' approach introduced by Sturmfels [15], can be employed to solve the linear system of {(Subscript[g, 1]), (Subscript[g, 2]), (Subscript[g, 3])}. 
Depending on which variable one wants, the original system is rewritten such that this particular variable is hidden (i.e. is treated as a constant). If our interest is to solve (Subscript[x, 1])= g((Subscript[x, 4])), the equations are first homogenized using a new variable (Subscript[x, 5]) and consider the variables  (Subscript[x, 1]) and (Subscript[x, 4]) as well as the constant part as parameters,

f1 = ((Subscript[a, 0, 3]) x1 + (Subscript[d, 3, 0]) x4 + (Subscript[e, 0, 3])) x5 + 
          (Subscript[b, 0, 3]) x2 + (Subscript[c, 0, 3]) x3;

f2 = ((Subscript[a, 1, 3]) x1  +  (Subscript[d, 3, 1]) x4 + (Subscript[e, 1, 3])) x5 + 
          (Subscript[b, 1, 3]) x2 +  (Subscript[c, 1, 3]) x3;

f3 = ((Subscript[a, 2, 3]) x1 + (Subscript[d, 3, 2]) x4 + (Subscript[e, 2, 3])) x5 + 
          (Subscript[b, 2, 3]) x2 + (Subscript[c, 2, 3]) x3;

The  Jacobian determinant becomes,

Jx1=Det[(({{(∂_x2) f1, (∂_x3) 
                    f1, (∂_x5) f1}, {(∂_x2) f2, (∂_x3) 
                    f2, (∂_x5) f2}, {(∂_x2) f3, (∂_x3) 
                    f3, (∂_x5) f3}}))];

Then the solution for (Subscript[x, 1]) as function of (Subscript[x, 4]) is,

solx1 = Solve[Jx1 == 0, x1]

Similarly, the homogenized system for  (Subscript[x, 2]) = g((Subscript[x, 4]))

f4 = ((Subscript[b, 0, 3]) x2 + (Subscript[d, 3, 0]) x4 + (Subscript[e, 0, 3])) x5 + 
          (Subscript[a, 0, 3]) x1 + (Subscript[c, 0, 3]) x3;

f5 = ((Subscript[b, 1, 3]) x2 +  (Subscript[d, 3, 1]) x4 + (Subscript[e, 1, 3])) x5 + 
          (Subscript[a, 1, 3]) x1  +  (Subscript[c, 1, 3]) x3;

f6 = ((Subscript[b, 2, 3]) x2 + (Subscript[d, 3, 2]) x4 + (Subscript[e, 2, 3])) x5 + 
          (Subscript[a, 2, 3]) x1 + (Subscript[c, 2, 3]) x3;

then

Jx2=Det[(({{(∂_x1) f4, (∂_x3) 
                    f4, (∂_x5) f4}, {(∂_x1) f5, (∂_x3) 
                    f5, (∂_x5) f5}, {(∂_x1) f6, (∂_x3) 
                    f6, (∂_x5) f6}}))];

and its solution is,

solx2 = Solve[Jx2 == 0, x2]

Finally (Subscript[x, 3]) =g((Subscript[x, 4])) leads to

f7 = ((Subscript[c, 0, 3]) x3 + (Subscript[d, 3, 0]) x4 + (Subscript[e, 0, 3])) x5 + 
          (Subscript[a, 0, 3]) x1 + (Subscript[b, 0, 3]) x2;

f8 = ( (Subscript[c, 1, 3]) x3 +  (Subscript[d, 3, 1]) x4 + (Subscript[e, 1, 3])) x5 + 
          (Subscript[a, 1, 3]) x1  + (Subscript[b, 1, 3]) x2;

f9 = ((Subscript[c, 2, 3]) x3 + (Subscript[d, 3, 2]) x4 + (Subscript[e, 2, 3])) x5 + 
          (Subscript[a, 2, 3]) x1 + (Subscript[b, 2, 3]) x2;

Jx3=Det[(({{(∂_x1) f7, (∂_x2) 
                    f7, (∂_x5) f7}, {(∂_x1) f8, (∂_x2) 
                    f8, (∂_x5) f8}, {(∂_x1) f9, (∂_x2) 
                    f9, (∂_x5) f9}}))];

solx3 = Solve[Jx3 == 0, x3]

On substituting the obtained expressions of  (Subscript[x, 1]) =g((Subscript[x, 4])), (Subscript[x, 2]) =g((Subscript[x, 4])) and (Subscript[x, 3]) =g((Subscript[x, 4])) in the fourth equation  (Subscript[e, 4]) ((Subscript[x, 1]), (Subscript[x, 2]), (Subscript[x, 3]), (Subscript[x, 4])),

G = 
        (Subscript[e, 4]) /. {solx1[[1, 1]], solx2[[1, 1]], 
            solx3[[1, 1]]};

This is a quadratic equation for  (Subscript[x, 4])

Exponent[G, x4, List]

The coefficients of this equation,  (Subscript[h, 2]) ((Subscript[x, 4])^2) + (Subscript[h, 1]) (Subscript[x, 4])  + (Subscript[h, 0]) = 0 are quite long expressions, therefore here we do not display them.

h2 = Coefficient[G, (x4^2)] // Simplify;

h1 = Coefficient[G, x4] // Simplify;

h0 = Simplify[G - (h2 (x4^2) + h1 x4)];

Indeed, this is a quadratic equation,

D[h0, x4]

The actual numeric solution of original system means the evaluation of  these coefficients with the numerical data,

h2c = h2 /. coeffsn /. data // Timing

h1c = h1 /. coeffsn /. data // Timing

h0c = h0 /. coeffsn /. data // Timing

and then solving the quadratic equation for variable  (Subscript[x, 4]),

solx4 = 
      Solve[h2c[[2]] x4^2 + h1c[[2]] x4 + h0c[[2]] == 0, 
          x4] // Timing

The two solutions are

{x41, x42} = x4 /. solx4[[2]]

Substituting these values in (Subscript[x, 1])= g((Subscript[x, 4])), we get (Subscript[x, 1])

{X1, X2} = 
      Map[x1 /. solx1[[1, 1]] /. coeffsn /. data /. 
            x4 -> # &, {x41, x42}]

Similarly, the values for (Subscript[x, 2]) are

{Y1, Y2} = 
      Map[x2 /. solx2[[1, 1]] /. coeffsn /. data /. 
            x4 -> # &, {x41, x42}]

and for  (Subscript[x, 3])

{Z1, Z2} = 
      Map[x3 /. solx3[[1, 1]] /. coeffsn /. data /. 
            x4 -> # &, {x41, x42}]

We can select the proper solution according to their norms,

Norm[{X1, Y1, Z1}]

Norm[{X2, Y2, Z2}]

Consequently, the first solution is correct,

SetPrecision[{X1, Y1, Z1}, 10]

It can be seen, that the direct numerical solution and the symbolic solution are the same for ten decimal digits (mm), while the direct numerical solution takes about two times more computing time, than  the symbolic solution.

In order to get the symbolic expression for the coefficients of the quadratic equation, other methods also can be used.

3.2.2 Dixon resultant

The Dixon resultant, Dixon [16] and its generalization for multivariate case Kapur et al. [17] are  elimination techniques similar to the Sylvester resultant for polynomial systems, Bernstein [18]. This function is not a built-in part of Mathematica but the package utilizes the algorithm of Kapur, Saxena and Yang, written by Nakos and Williams [19] is avaiable from Wolfram Mathsource, on the web at: http://library.wolfram.com/infocenter/Articles/2597/.

<< Resultant`Dixon`

Now we can solve the original system, eliminating the variables  (Subscript[x, 2]) and (Subscript[x, 3]) to get  (Subscript[x, 1])= g((Subscript[x, 4])),

drx1 = 
      DixonResultant[{(Subscript[g, 1]), (Subscript[g, 2]), (Subscript[g, 3])}, {x2, x3}, {u2, u3}]

This is a linear expression containing only  (Subscript[x, 1]) and (Subscript[x, 4]),

Exponent[drx1, {x1, x2, x3, x4}]

Then the solution  for (Subscript[x, 1])= g((Subscript[x, 4]))

solx1 = Solve[drx1 == 0, x1]

Similarly, for the two additional variables, (Subscript[x, 2])= g((Subscript[x, 4])) and (Subscript[x, 3])= g((Subscript[x, 4])),

drx2 = 
      DixonResultant[{(Subscript[g, 1]), (Subscript[g, 2]), (Subscript[g, 3])}, {x1, x3}, {u1, u3}]

Exponent[drx2, {x1, x2, x3, x4}]

solx2 = Solve[drx2 == 0, x2]

and

drx3 = 
      DixonResultant[{(Subscript[g, 1]), (Subscript[g, 2]), (Subscript[g, 3])}, {x1, x2}, {u1, u2}]

Exponent[drx3, {x1, x2, x3, x4}]

solx3 = Solve[drx3 == 0, x3]

After substitution, we have again a quadratic equation for (Subscript[x, 4]),

G = 
        (Subscript[e, 4]) /. {solx1[[1, 1]], solx2[[1, 1]], 
            solx3[[1, 1]]};

Exponent[G, x4, List]

The coefficents of the quadratic equation are,

h2d = Coefficient[G, (x4^2)];

h1d = Coefficient[G, x4];

h0d = Simplify[G - (h2d (x4^2) + h1d x4)];

and

D[h0, x4]

The coefficients provided by the Sturmfels' method and the Dixon resultant are the same,

{h2 - h2d, h1 - h1d, h0 - h0d} // Simplify

3.2.3 Groebner basis approach

Probably, the most effective elimination technique is the lexicographic Groebner basis utilizes the Gauss elimination technique for multivariate polynomial systems proposed by Buchberger [20]. This method is implemented in Mathematica and Maple, too.

gb = 
      GroebnerBasis[{(Subscript[g, 1]), (Subscript[g, 2]), (Subscript[g, 3])}, {x1, x2, x3, x4}] // 
        Simplify

This basis consists of seven polynoms,

Length[gb]

When we compute the exponents of the variables of the different polynomials,

Table[Exponent[gb[[i]], {x1, x2, x3, x4}], {i, 1, 7}]

it is easy to see, that from the first one, one can get (Subscript[x, 3])= g((Subscript[x, 4])), then employing this,  the second one gives (Subscript[x, 2])= g((Subscript[x, 4])) and the fifth one provides (Subscript[x, 1])= g((Subscript[x, 4])). Let us compute (Subscript[x, 3]) as function of (Subscript[x, 4]),

solx3 = Solve[gb[[1]] == 0, x3] // Simplify

Then  (Subscript[x, 2])= g((Subscript[x, 4]))

solx2 = 
      Solve[gb[[2]] == 0, x2] /. solx3[[1]] // 
        Simplify

Finally, (Subscript[x, 1])= g((Subscript[x, 4]))

solx1 = Solve[gb[[5]] == 0, x1] /. solx2[[1]] /. 
          solx3[[1]] // Simplify

After substitution,  we get quadratic equation again,

G = 
        (Subscript[e, 4]) /. {solx1[[1, 1]], solx2[[1, 1]], 
            solx3[[1, 1]]};

Exponent[G, x4, List]

Then, we get the coefficients of the quadratic equation in the usual way,

h2gL = Coefficient[G, (x4^2)];

h1gL = Coefficient[G, x4];

h0gL = Simplify[G - (h2gL (x4^2) + h1gL x4)];

D[h0gL, x4]

This solution is correct indeed,

{h2 - h2gL, h1 - h1gL, h0 - h0gL} // Simplify

3.2.4 Reduced Groebner basis approach

The built-in function GroebnerBasis can be also used to eliminate directly prespecified variables  form the basis, Trott [21]. In this case, one may employ this function in similar way than the Dixon method does. It goes without saying, that now the function should be called three times with different arguments.

First, again we want (Subscript[x, 1])= g((Subscript[x, 4])), therefore variables (Subscript[x, 2]) and (Subscript[x, 3]) should be eliminated from the Groebner basis,

{gbx1} = 
        GroebnerBasis[{(Subscript[g, 1]), (Subscript[g, 2]), (Subscript[g, 3])}, {x1, x4}, {x2, x3}];

Now, the basis contains only one equation,

Length[{gbx1}]

in which only (Subscript[x, 1]) and (Subscript[x, 4]) can be found,

Exponent[gbx1, {x1, x2, x3, x4}]

Therefore (Subscript[x, 1])  =  
      g ((Subscript[x, 4])) can be computed directly,

solx1 = Solve[gbx1 == 0, x1] // Simplify

Similarly, in the other cases

{gbx2} = 
        GroebnerBasis[{(Subscript[g, 1]), (Subscript[g, 2]), (Subscript[g, 3])}, {x2, x4}, {x1, x3}];

Exponent[gbx2, {x1, x2, x3, x4}]

solx2 = Solve[gbx2 == 0, x2] // Simplify

and

{gbx3} = 
        GroebnerBasis[{(Subscript[g, 1]), (Subscript[g, 2]), (Subscript[g, 3])}, {x3, x4}, {x1, x2}];

Exponent[gbx3, {x1, x2, x3, x4}]

solx3 = Solve[gbx3 == 0, x3] // Simplify

After substituting them, we get

G = 
        (Subscript[e, 4]) /. {solx1[[1, 1]], solx2[[1, 1]], 
            solx3[[1, 1]]};

Exponent[G, x4, List]

Then the coefficients of the quadratic equation are,

h2gR = Coefficient[G, (x4^2)];

h1gR = Coefficient[G, x4];

h0gR = Simplify[G - (h2gR (x4^2) + h1gR x4)];

D[h0gR, x4]

We have again, the same result,

{h2 - h2gR, h1 - h1gR, h0 - h0gR} // Simplify

The values of  the input data are seldom known exactly because of measurement errors. Consequently, the following question can arise: How does the relative error in the solution depend on the error of the input data? The sensitivity of a solution can be estimated on the bases of the linearized form of the system of the polynomial equations ({e1, e2, e3, e4}, around the solution. The bound of the relative error approximately proportional with the condition number of the matrix of this linearized system as well as with working precision of the computation, see Skeel and Keiper [27]. 
In practice, in order to reduce measurement error, the signals of more than four satellites have been received.

4.  Solutions of  the N-points problem

Now, we extend the concept to the case  where more than four GPS satellites are in view as it is usually in practice. In that case, there are more independent equations than unknown variables.  For illustration, we have data of six satellites with the following numerical values,

datan = {(Subscript[a, 0]) -> 14177553.47, (Subscript[a, 1]) -> 15097199.81, 
          (Subscript[a, 2]) -> 23460342.33, (Subscript[a, 3]) -> -8206488.95, 
          (Subscript[a, 4]) -> 1399988.07, 
          (Subscript[a, 5]) -> 
            6995655.48, 
(Subscript[b, 0]) -> -18814768.09, 
          (Subscript[b, 1]) -> -4636088.67, 
          (Subscript[b, 2]) -> -9433518.58, 
(Subscript[b, 3]) -> -18217989.14, (Subscript[b, 4]) -> -17563734.90, 
          (Subscript[b, 5]) -> -23537808.26, 
(Subscript[c, 0]) -> 
            12243866.38, (Subscript[c, 1]) -> 21326706.55, 
          (Subscript[c, 2]) -> 8174941.25, 
(Subscript[c, 3]) -> 
            17605231.99, (Subscript[c, 4]) -> 19705591.18, 
          (Subscript[c, 5]) -> -9927906.48,  
(Subscript[d, 0]) -> 
            21119278.32, (Subscript[d, 1]) -> 22527064.18, 
          (Subscript[d, 2]) -> 23674159.88, 
(Subscript[d, 3]) -> 
            20951647.38, (Subscript[d, 4]) -> 20155401.42, 
          (Subscript[d, 5]) -> 24222110.91};

First, as direct numerical method, the minimization of the residiums of the equation, the least square approach will be presented.

4.1  Direct numerical solution

The general form of the equation for the i-th satellite is

e = (((x1 - (Subscript[a, i])))^2) + (((x2 - (Subscript[b, i])))^2) + (((x3 - (Subscript[c, i])))^2) - (((x4 - (Subscript[d, i])))^2);

The objective to be minimized is the sum of the residium of the equations

f = Apply[Plus, Table[(e^2) /. datan, {i, 0, 5}]] // Simplify

To find the global minimum, the built-in function NMinimize is more promising than FindMinimum designed to search local minimum,

solN = NMinimize[f, {x1, x2, x3, x4}] // Timing

or

SetPrecision[solN[[2, 2]], 10]

However, if we employ the norm of the distance instead of the residium of the equations, namely

en = ((((Subscript[d, i]) - ((Sqrt[(((x1 - (Subscript[a, i])))^2) + (((x2 - (Subscript[b, i])))^2) + (((x3 - (Subscript[c, i])))^2)]) - x4)))^2);

and then the objective is,

fn = (Sqrt[Apply])[Plus, Table[en /. datan, {i, 0, 5}]]

The optimum will be somewhat different,

solNn = NMinimize[fn, {x1, x2, x3, x4}] // Timing

or

SetPrecision[solNn[[2, 2]], 10]

4.2  Indirect semi-symbolical solution

This minimization problem can be transformed into a set of nonlinear equations. The symbolic form of the objective to be minimized is,

f = Apply[Plus, Table[(e^2), {i, 0, 5}]]

The necessary conditions for the existence of the extreme value  are the zero partial derivatives, namely

g1 = D[f, x1];

g2 = D[f, x2];

g3 = D[f, x3];

g4 = D[f, x4];

In order to solve this nonlinear polynomial  system,  the Groebner basis approach can be used. However to compute the Groebner basis in fully symbolic form, would take a very long computation time, and the expressions would be very long providing no extra advantages. Therefore, here we compute the analytical form of the reduced Groebner basis, namely for (Subscript[x, 1])

gbx1 = 
      GroebnerBasis[{g1, g2, g3, g4} /. datan, {x1}, {x2, x3, x4}]

which means, that no Groebner basis exists on this low precision coefficent domain, consequently higher precision is needed. The rationalization of the data ensures infinite precision,

datanR = 
      Map[#[[1]] -> Rationalize[#[[2]]] &, datan]

However, we can considerably reduce computation time by employing only high-precision numbers instead of rational numbers. Now the reduced  Groebner basis for (Subscript[x, 1]), on this domain is

gbx1 = GroebnerBasis[
            N[{g1, g2, g3, g4} /. datanR, 100], {x1}, {x2, x3, x4}, 
            CoefficientDomain -> InexactNumbers] // N // 
        Timing

Solving this polynomial equation,

Solx1 = NSolve[gbx1[[2]] == 0, x1]; // Timing

Selecting the positive real solutions, we get

x1s = 
      Select[Map[x1 /. # &, Flatten[Solx1]], Im[#] == 0 &]

Similarly,  for variables  (Subscript[x, 2]),  (Subscript[x, 3]) and  (Subscript[x, 4])

gbx2 = 
        GroebnerBasis[N[{g1, g2, g3, g4} /. datanR, 100], {x2}, {x1, x3, x4}, 
          CoefficientDomain -> InexactNumbers];

Solx2 = NSolve[gbx2 == 0, x2];

x2s = 
      Select[Map[x2 /. # &, Flatten[Solx2]], Im[#] == 0 &] // 
        N

gbx3 = 
        GroebnerBasis[N[{g1, g2, g3, g4} /. datanR, 100], {x3}, {x1, x2, x4}, 
          CoefficientDomain -> InexactNumbers];

Solx3 = NSolve[gbx3 == 0, x3];

x3s = 
      Select[Map[x3 /. # &, Flatten[Solx3]], Im[#] == 0 &] // 
        N

gbx4 = 
        GroebnerBasis[N[{g1, g2, g3, g4} /. datanR, 100], {x4}, {x1, x2, x3}, 
          CoefficientDomain -> InexactNumbers];

Solx4 = NSolve[gbx4 == 0, x4];

x4s = 
      Select[Map[x4 /. # &, Flatten[Solx4]], Im[#] == 0 &] // 
        N

In order to decide, which solution is acceptable, the norm of distance vector from the geometrical coordinates, {(Subscript[x, 1]),(Subscript[x, 2]), (Subscript[x, 3])} can be computed,

x123 = Transpose[{x1s, x2s, x3s}];

Map[Norm[#] &, x123]

So, we can see that the second solution is correct,

SetPrecision[x123[[2]], 10]

The fourth variable, (Subscript[x, 4]) can be computed from one of the equations of the system, let say from the first one,

g1 /. datanR /. {x1 -> x123[[2, 1]], x2 -> x123[[2, 2]], 
          x3 -> x123[[2, 3]]} // Simplify

A correct value is, which can be found in the Groebner solution, too (see above x4s), therefore the acceptable value is (Subscript[x, 4]) = 13.4526.
So we have got the same solution as the result of the direct numerical minimization.

4.3  Gauss-Jacobi combinatorical solution

The Gauss-Jacobi combinatorical solution, Jacobi [22] can be employed originally for linear regression problem.  This means, if we have more independent equations, m than variables, n, so m > n, then the solution, in least-squares sense, can be achieved by solving the  (({{m}, {n}})) combinatorical subsets of n equations, and then weighting these solutions properly, Wellisch [23] and  Hornoch [24]. This method can be extended to nonlinear case using linearization. A key point of this extention, that we can solve nonlinear subsets without a priori information about the solutions, which would be necessary for the traditional linearization. The determination of weights through linearization takes place after the solutions of these nonlinear subsets are known, see for example Závoti [25].

First, the subsets should be determined. In our case

n = 4; m = 6;

The number of the subsets

Binomial[m, n]

These subsets  are,

qs = 
      Partition[Map[# - 1 &, Flatten[Subsets[Range[m], {n}]]], 4]

The values of  the indices start from zero in correspondence of the indices of the coefficients of the equations. Now, we shall utilize the symbolic solution of  the GPS 4-points problem, namely the expressions of the coefficients of the quadratic equation ((Subscript[h, 2]),(Subscript[h, 1]),(Subscript[h, 0])). Therefore, we construct a new data list, datap similar to datan,which assignes the proper values to the  coefficients of the equations of the 15 subsets,

datap = 
        Table[Map[
            Select[datan,  
                  MemberQ[
                      qs[[i]], #[[1, 2]]] &] /. {#[[1]] -> 
                    0, #[[2]] -> 1, #[[3]] -> 
                    2, #[[4]] -> 3} &, {qs[[i]]}], {i, 1, 
            15}];

For example, the fourth subset is indexed as

qs[[4]]

and it has the proper data assignments,

datap[[4]]

Now, we can employ the symbolic expressions of the coefficients of the quadratic equation for (Subscript[x, 4]), ((Subscript[h, 2]), (Subscript[h, 1]), (Subscript[h, 0])), which were developed for the GPS 4-points problem, see 3.2.1. These coefficients can be evaluated for all of the 15 combinatorical subsets, ((Subscript[H, 2]), (Subscript[H, 1]), (Subscript[H, 0])),

H2 = Map[(h2 /. coeffsn /. Flatten[#]) &, datap]; // 
      Timing

H1 = Map[(h1 /. coeffsn /. Flatten[#]) &, datap]; // 
      Timing

H0 = Map[(h0 /. coeffsn /. Flatten[#]) &, datap]; // 
      Timing

It is useful to display these coeffients,

H210 = Transpose[{H2, H1, H0}];

H210c = Map[# /. datan &, H210]; 
    TableForm[SetPrecision[H210c, 15]]

This table indicates that the 10th  combination has a poor geometry which fact can be also detected by computing its PDOP (Position  Dilution of Precision), see Hofman et al. [26].
Then the 15 quadratic equations can be solved for (Subscript[x, 4]),

X4 = 
          Map[x4 /. Solve[#[[1]] x4^2 + #[[2]] x4 + #[[3]] == 0, x4][[1, 1]] &, H210c]; // Timing

These values of (Subscript[x, 4]) can be substituted into the symbolic relations (Subscript[x, 1]) = g((Subscript[x, 4])), (Subscript[x, 2]) = g((Subscript[x, 4])) and (Subscript[x, 3]) = g((Subscript[x, 4])) developed for GPS 4 - points problem, see 3.2.4,

X1 = 
        MapThread[(x1 /. solx1[[1, 1]] /. coeffsn /. 
                  Flatten[#1] /. x4 -> #2) &, {datap, 
            X4}];

X2 = 
        MapThread[(x2 /. solx2[[1, 1]] /. coeffsn /. 
                  Flatten[#1] /. x4 -> #2) &, {datap, 
            X4}];

X3 = 
        MapThread[(x3 /. solx3[[1, 1]] /. coeffsn /. 
                  Flatten[#1] /. x4 -> #2) &, {datap, 
            X4}];

Let us display these solutions for  ((Subscript[x, 1]), (Subscript[x, 2]), (Subscript[x, 3]), (Subscript[x, 4])),

X = Transpose[{X1, X2, X3, X4}]; 
    TableForm[SetPrecision[X, 11]]

In order to compute the weights of these solutions, one has to compute the square of the 15 Jacobi determinants. Each has size of  4 × 4, because of four equations and four variables.  Starting with the general form of the i-th equation,

e

The partial derivatives are,

de = {D[e, x1], D[e, x2], D[e, x3], D[e, x4]}

The numerical values of these partial derivatives will be computed at the corresponding combinatorical solutions, see the table above (variable X). Therefore, the weights, (Subscript[π, j]),  the square of the 15 Jacobi determinants are,

πs = 
          Table[Map[(((Det[{(de /. i -> #[[1]]), (de /. i -> #[[2]]), (de /. i -> #[[3]]), (de /. i -> #[[4]])} /. datan]))^2) &, qs][[j]] /. {x1 -> 
                  X[[j, 1]], x2 -> X[[j, 2]], 
                x3 -> X[[j, 3]], x4 -> X[[j, 4]]}, {j, 1, 
              15}]; // Timing

The sum of these weights are sπ,

sπ = (∑_(j = 1)^15) (Subscript[π, j])

sπs = Apply[Plus, πs]

Then the weighted solution of the variable (Subscript[x, i]) is

(Subscript[x, i])=(1/s)π (∑_(j = 1)^15)(((Subscript[x, i]))^((j)))

{X1s, X2s, X3s, X4s} = Map[πs . # &, {X1, X2, X3, X4}]/sπs; 
    SetPrecision[{X1s, X2s, X3s, X4s}, 10]

which is  very close to the direct numerical solution, which was

Conclusions

In Table I. we summarized the computation times of the different methods. The computations were carried out with Mathematica 5.2 on HP xw 4100 workstation employing Pentium 3 MHz processor and 1 MB RAM.

Table I. Computation times in seconds

({{Method, 4 - points, N - 
              points (N  =  6)}, {Direct numerial, 0.078, 0.094}, {Symbolical, 0.031, -}, {Semi - symbolical, -, 0.188}, {Gauss-Jacobi , -, 0.718}})

In case of 4-points problem the symbolical solution proved to be more than two times faster than the direct numerical one. For N-points problem (N = 6) the direct numerical is about two times faster than the semi - symbolical solution both of them are based on numerical minimization, while the Gauss - Jacobi solution needs nearly eight times more computation time than the direct numerical solution. However, the latest one can provide very important information about poor geometry or outliers in observations, Awange and Grafarend [2].

Acknowledgements

The author gratefully acknowledge the suggestions of the referees for improving this paper.

References

[1] Haneberg, W.C.:Computational Geosciences with Mathematics, Springer , Berlin, (2004)

[2] Awange J.L. and Grafarend E.W.: Solving Algebraic Computational Problems in Geodesy and Geoinformatics, Springer, Berlin, (2005)

[3] Freeman, J.A.: Simulating neural networks with Mathematica, Addison-Wesley, (1994)

[4] Bart M. ter Haar Romeny: Front-End Vision and Multi-Scale Image Analysis, Kluwer, Dordrecht,  (2003)

[5] Helton, J.W and Merino O.: Classical Control using (H^∞) Methods, SIAM,  (1998)

[6] Bellomo, N., Preziosi, L. and Romano, A.: Mechanics and Dynamical Systems with Mathematica, Birkhäuser, (2000)

[7] Bancroft, S.: An algebraic solution of the GPS equations, IEEE Transaction on Aerospace and Electronic Systems AES-21, (1985), pp. 56-59.

[8] Grafarend EW. and Schaffrin, B.: Ausglechungsrechnung in Linearen Modellen, B.I. Wissenschaftsverlag, Mannheim,  (1993) 

[9] Singer P., Ströbel D., Hördt R., Bahndorf J., Linkwitz K.: Direkte Lösung des räumlichen Bogenschnitts. Zeitschrift für Vermessungswesen, 124, (1993), S. 295 - 297.

[10] Kleusberg A.: Die direkte Lösung des räumlichen Hyperbelschnitts. Zeitschrift für Vermessungswesen, 119,  (1994), S. 188 - 192.

[11] Kleusberg  A.: Analytical GPS navigation solution. In Grafarend EW., Krumm FW. and Schwarze VS. (eds) Geodesy -the Challenge of the 3rd Millennium, Springer, Heidelberg, (2003), pp. 93 - 96.

[12] Lichtenegger H.: Eine direkte Lösung des räumlichen Bogenschnitts., Österreichische Zeitschrift für Vermessung und Geoinformation 83,  (1995), S. 224 - 226.

[13] Awange J.L.: Gröbner bases, multipolynomial resultants and Gauss-Jacobi combinatorial algorithms - adjustment of nonlinear GPS/LPS observations, Dissertation (D93), Geodätisches Institute der Universität Stuttgart, (2002).

[14] Grafarend EW. and Shan J.: Closed - form solution of the nonlinear  pseudoranging equations (GPS), Artificial Satellites, Planetary Geodesy 31,  (1996), pp.133 - 147.

[15] Sturmfels B.: Multigraded resultant of  Sylvester type. Journal of Algebra 163, (1994), pp.115-127.

[16] Dixon A.L.: The eliminant of three quantics in two independent variables. Proc. London Math.Soc. 6,  (1908),  pp. 468-478.

[17] Kapur D, Saxena T. and Yang L. (1994): Algebraic and geometric reasoning using Dixon resultants. In: ACM ISSAC 94. Oxford, England, July 1994, pp. 99-107.

[18] Bernstein, D.S.: Matrix Mathematics, Theory, Facts, and Formulas with Application to Linear Systems Theory, Princeton University Press, Princeton and Oxford, (2005), p.142.

[19] Nakos G. and Williams R.: Elimination with The Dixon resultant. Mathematica in Education and Research, Vol.6. No.2,  (1997), pp. 11 - 21.

[20] Buchber B. Gröbner bases: An Algorithmic method in Polynomial Ideal theory, Multidimensional Systems Theory, N.K. Bose, ed., D.Reidel Publ.Co. (1985)

[21] Trott M.: Mathematica Guide Book for Symbolic Computation, Springer, New York,  (2006)

[22] Jacobi CGI.: Deformatione et proprietatibus determinantum, Crelle's Journal für die reine und angewandte Mathematik, Bd. 22,  (1841)

[23]  Wellisch S.: Theorie und Praxis der Ausgleichungsrechnung. Bd.2, (1910)
 
[24] Tárczy - Hornoch A.: Über die Zurückführung der Methode der kleinsten Quadrate  auf das Prinzipe des arithmetischen Mittels. Österreichische Zeitschrift für Vermessungswesen 38, (1950),  13 - 18.
  
[25] Závoti J.: The closed form solution of  7 parameter 3D transformation. Publications in Geomatics, Sopron, Hungary, (2005), pp. 53 -60.

[26] Hofman-Wellenhof B, Lichtenegger H, Collins J.: Global Positioning System: Theory and practice, 5th Ed. Springer, Wien, (2001).

[27] Skeel R.D. and Keiper J.B. Elementary Numerical Computing With Mathematica, McGraw-Hill, (1999)