CODEWITHSUNDEEP
pythonmathscipymathematical-optimizationscipy-optimize
10GeV
10GeV

Reputation: 475

Solving this rectangular, nonlinear system with SciPy

Background.

I'm attempting to write a implementation of this answer over on Math SE. You may find the following background to be useful.

Problem

I have an experimental setup consisting of three (3) receivers, with known locations [xi, yi, zi], and a transmitter with unknown location [x,y,z] emitting a signal at known velocity v. This signal arrives at the receivers at known times ti. The time of emission, t, is unknown.

I wish to find the angle of arrival (i.e. the transmitter's polar coordinates theta and phi), given only this information.

Solution

It is not possible to locate the transmitter exactly with only three (3) receivers, except in a handful of unique cases (there are several great answers across Math SE explaining why this is the case). In general, at least four (and, in practice, >>4) receivers are required to uniquely determine the rectangular coordinates of the transmitter.

The direction to the transmitter, however, may be "reliably" estimated. Letting vi be the vector representing the location of receiver i, ti being the time of signal arrival at receiver i, and n be the vector representing the unit vector pointing in the (approximate) direction of the transmitter, we obtain the following equations:

<n, vj - vi> = v(ti - tj)

(where < > denotes the scalar product)

...for all pairs of indices i,j. Together with |n| = 1, the system has 2 solutions in general, symmetric by reflection in the plane through vi/vj/vk. We may then determine phi and theta by simply writing n in polar coordinates.

Implementation.

I've attempted to write a implementation of the above solution, using 's fsolve.

from dataclasses import dataclass
import scipy.optimize
import random
import math

c = 299792

@dataclass
class Vertexer:

        roc: list

        def fun(self, var, dat):
                (x,y,z) = var

                eqn_0 = (x * (self.roc[0][0] - self.roc[1][0])) + (y * (self.roc[0][1] - self.roc[1][1])) + (z * (self.roc[0][2] - self.roc[1][2])) - c * (dat[1] - dat[0])
                eqn_1 = (x * (self.roc[0][0] - self.roc[2][0])) + (y * (self.roc[0][1] - self.roc[2][1])) + (z * (self.roc[0][2] - self.roc[2][2])) - c * (dat[2] - dat[0])
                eqn_2 = (x * (self.roc[1][0] - self.roc[2][0])) + (y * (self.roc[1][1] - self.roc[2][1])) + (z * (self.roc[1][2] - self.roc[2][2])) - c * (dat[2] - dat[1])

                norm = math.sqrt(x**2 + y**2 + z**2) - 1

                return [eqn_0, eqn_1, eqn_2, norm]

        def find(self, dat):
                result = scipy.optimize.fsolve(self.fun, (0,0,0), args=dat)
                print('Solution ', result)

# Crude code to simulate a source, receivers at random locations

x0 = random.randrange(0,50); y0 = random.randrange(0,50); z0 = random.randrange(0,50)

x1 = random.randrange(0,50); x2 = random.randrange(0,50); x3 = random.randrange(0,50);
y1 = random.randrange(0,50); y2 = random.randrange(0,50); y3 = random.randrange(0,50);
z1 = random.randrange(0,50); z2 = random.randrange(0,50); z3 = random.randrange(0,50);

t1 = math.sqrt((x0-x1)**2 + (y0-y1)**2 + (z0-z1)**2)/c
t2 = math.sqrt((x0-x2)**2 + (y0-y2)**2 + (z0-z2)**2)/c
t3 = math.sqrt((x0-x3)**2 + (y0-y3)**2 + (z0-z3)**2)/c

print('Actual coordinates ', x0,y0,z0)

myVertexer = Vertexer([[x1,y1,z1], [x2,y2,z2], [x3,y3,z3]])
myVertexer.find([t1,t2,t3])

Unfortunately, I have far more experience solving such problems in C/C++ using GSL, and have limited experience working with and the like. I'm getting the error:

TypeError: fsolve: there is a mismatch between the input and output shape of the 'func' argument 'fun'.Shape should be (3,) but it is (4,).

...which seems to suggest that fsolve expects a square system.

How may I solve this rectangular system? I can't seem to find anything useful in the docs.

If necessary, I'm open to using other (Python) libraries.

Upvotes: 0

Views: 176

Answers (1)

joni
joni

Reputation: 7157

As you already mentioned, fsolve expects a system with N variables and N equations, i.e. it finds a root of the function F: R^N -> R^N. Since you have four equations, you simply need to add a fourth variable. Note also that fsolve is a legacy function, and it's recommended to use root instead. Last but not least, note that sqrt(x^2+y^2+z^2) = 1 is equivalent to x^2+y^2+z^2=1 and that the latter is much less susceptible to rounding errors caused by the finite differences when approximating the jacobian of F.

Long story short, your class should look like this:

from scipy.optimize import root

@dataclass
class Vertexer:
        roc: list
        def fun(self, var, dat):
                x,y,z, *_ = var
                eqn_0 = (x * (self.roc[0][0] - self.roc[1][0])) + (y * (self.roc[0][1] - self.roc[1][1])) + (z * (self.roc[0][2] - self.roc[1][2])) - c * (dat[1] - dat[0])
                eqn_1 = (x * (self.roc[0][0] - self.roc[2][0])) + (y * (self.roc[0][1] - self.roc[2][1])) + (z * (self.roc[0][2] - self.roc[2][2])) - c * (dat[2] - dat[0])
                eqn_2 = (x * (self.roc[1][0] - self.roc[2][0])) + (y * (self.roc[1][1] - self.roc[2][1])) + (z * (self.roc[1][2] - self.roc[2][2])) - c * (dat[2] - dat[1])
                norm = x**2 + y**2 + z**2 - 1
                return [eqn_0, eqn_1, eqn_2, norm]

        def find(self, dat):
                result = root(self.fun, (0,0,0,0), args=dat)
                if result.success:
                        print('Solution ', result.x[:3])

Upvotes: 2

Related Questions