Problem.
Consider the following problem: An experiment consists of three receivers, i,j and k, with known coordinates (xi,yi,zi), (xj,yi... in 3D space. A stationary transmitter, with unknown coordinates (x,y,z), emits a signal with known velocity v. The time of arrival of this signal at each receiver is recorded, giving ti,tj,tk. The signal emission time, t, is unknown.
Using only the coordinates of the receivers and the signal arrival times, I wish to determine the location of the transmitter.
In general, to solve for the coordinates of such a transmitter in N-dimensional space, N+1 receivers are required. Hence, in this case, a single, unique solution is unobtainable. A small number of finite solutions should be obtainable via numerical methods, however.
We may write the following system of equations to model the problem:
Sqrt[(x-xi)^2 + (y-yi)^2 + (z-zi)^2] + s(tj-ti) = Sqrt[(x-xj)^2 + (y-yj)^2 + (z-zj)^2]
Sqrt[(x-xj)^2 + (y-yj)^2 + (z-zj)^2] + s(tk-tj) = Sqrt[(x-xk)^2 + (y-yk)^2 + (z-zk)^2]
Sqrt[(x-xi)^2 + (y-yi)^2 + (z-zi)^2] + s(tk-ti) = Sqrt[(x-xk)^2 + (y-yk)^2 + (z-zk)^2]
Each equation gives a hyperboloid. Under ideal conditions, these three hyperboloids will intersect at precisely two points— one being the "true" solution, and the other being a reflection of that solution about the plane defined by the three receivers. In practice, given sufficiently accurate measurements, numerical solvers should be able to approximate these points of intersection.
My goal is to determine both solutions. Although it is impossible to determine which is the "true" transmitter location, for my purposes this will be sufficient.
Implementation.
I wish to implement a solution in Python. I'm not terribly familiar with NumPy or SciPy, however I've done a fair bit of work in SymPy, so I began there.
SymPy offers a variety of solvers, most of which focus on obtaining solutions symbolically. Not surprisingly, solve() and the like failed to find a solution, even under simulated "ideal" conditions (picking a random point, calculating the time taken for a signal originating from that point to arrive at each receiver, and feeding this to the algorithm).
SymPy also offers a numerical solver, nsolve(). I gave this a go, using the approach given below, however (not surpisingly) I got the error ZeroDivisionError: Matrix is numerically singular.
f = sym.Eq(sym.sqrt((x - x_i)**2 + (y - y_i)**2 + (z - z_i)**2) - sym.sqrt((x - x_j)**2 + (y - y_j)**2 + (z - z_j)**2), D_ij)
g = sym.Eq(sym.sqrt((x - x_i)**2 + (y - y_i)**2 + (z - z_i)**2) - sym.sqrt((x - x_k)**2 + (y - y_k)**2 + (z - z_k)**2), D_ik)
h = sym.Eq(sym.sqrt((x - x_j)**2 + (y - y_j)**2 + (z - z_j)**2) - sym.sqrt((x - x_k)**2 + (y - y_k)**2 + (z - z_k)**2), D_jk)
print("Soln. ", sym.nsolve((f,g,h),(x,y,z), (1,1,1)))
As I understand it, nsolve() relies on "matrix" techniques. A singular matrix is one which may not be inverted (or, equivalently, which has determinant zero), hence SymPy is unable to solve the system in question.
My understanding of matrices and nonlinear system solving techniques is a bit lacking, however my understanding is that a singular matrix occurs when there are infinitely many solutions, no solutions, or more than one solution. Given that I know there to be exactly two solutions, I believe this is the issue.
What Python solvers are available can be used so solve a nonlinear system with multiple solutions? Or, alternatively, is there a way to modify this such that it is digestible by SymPy?
Browsing the Numpy and SciPy docs, it seems that their solvers are effectively identical to what SymPy offers.
The (crude) test code mentioned below:
from random import randrange
import math
# SIMPLE CODE FOR SIMULATING A SIGNAL
# Set range
N=10
P=100
# Pick nodes to be at random locations
x_1 = randrange(N); y_1 = randrange(N); z_1 = randrange(N)
x_2 = randrange(N); y_2 = randrange(N); z_2 = randrange(N)
x_3 = randrange(N); y_3 = randrange(N); z_3 = randrange(N)
# Pick source to be at random location
x = randrange(P); y = randrange(P); z = randrange(P)
# Set velocity
c = 299792 # km/ns
# Generate simulated source
t_1 = math.sqrt( (x - x_1)**2 + (y - y_1)**2 + (z - z_1)**2 ) / c
t_2 = math.sqrt( (x - x_2)**2 + (y - y_2)**2 + (z - z_2)**2 ) / c
t_3 = math.sqrt( (x - x_3)**2 + (y - y_3)**2 + (z - z_3)**2 ) / c
# Normalize times to remove information about 'true' emission time
earliest = min(t_1, t_2, t_3)
t_1 = t_1 - earliest; t_2 = t2 - earliest; t_3 = t_3 = earliest