用 SciPy 求解这个矩形非线性系统

Question

背景。

我正在尝试在 Math SE 上写一个 python implementation of this 答案。您可能会发现以下背景很有用。

问题

我有一个实验装置，由三 (3) 个接收器组成，位置已知 [xi, yi, zi]，发射器的位置未知 [x,y,z] 已知速度v的信号。该信号在已知时间 ti 到达接收器。发射时间 t 未知.

我希望找到到达角（即发射机的极坐标theta和phi），仅给出此信息。

解决方案

不可能仅用三 (3) 个接收器就可以准确定位发射器，少数特殊情况除外（[=78= 中有几个很好的答案） ]Math SE 解释了为什么会这样）。一般情况下，至少需要四个（实际上是>>4个）接收器来唯一确定发射器的直角坐标。

然而，

方向可以“可靠地”估计到发射器。设 vi 是表示接收器位置的向量 i，ti 是信号到达接收器的时间 i，n 是表示接收器位置的向量指向发射机（近似）方向的单位向量，我们得到以下等式：

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

（其中< >表示标量积）

...对于所有索引对 i、j。连同 |n| = 1，该系统通常有 2 个解决方案，通过 vi/vj/vk 在平面中通过反射对称。然后我们可以通过简单地在极坐标中写入n来确定phi和theta。

实施。

我试图写一个 python implementation of the above solution, using scipy 的 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])

不幸的是，我在 C/C++ 中使用 GSL 解决此类问题的经验要丰富得多，而在 scipy 等方面的经验有限。我收到错误：

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,).

...这似乎表明 fsolve 需要一个方形系统。

我该如何解决这个矩形系统？我似乎在 scipy 文档中找不到任何有用的东西。

如有必要，我愿意使用其他 (Python) 库。

Answer 1

正如您已经提到的，fsolve 需要一个具有 N 个变量和 N 个方程的系统，即它找到函数 F 的根：R^N -> R^N。由于您有四个方程式，因此您只需添加第四个变量即可。另请注意，fsolve 是遗留函数，建议改用 root。最后但同样重要的是，请注意 sqrt(x^2+y^2+z^2) = 1 等同于 x^2+y^2+z^2=1 并且后者在近似 F.

的雅可比矩阵时更不容易受到有限差分引起的舍入误差的影响

长话短说，您的 class 应该如下所示：

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])

用 SciPy 求解这个矩形非线性系统

Solving this rectangular, nonlinear system with SciPy

python

math

scipy

mathematical-optimization

scipy-optimize

背景。

实施。