使用 NumPy 最小化此错误函数
Minimizing this error function, using NumPy
背景
我一直在努力解决(众所周知的痛苦)到达时间差 (TDoA) 多点定位问题,在 3-dimensions
中并使用 4
节点。如果你不熟悉这个问题,它是确定一些信号源的坐标(X,Y,Z)
,给定n
个节点的坐标,信号到达每个节点的时间,以及速度信号 v
.
我的解决方法如下:
对于每个节点,我们写(X-x_i)**2 + (Y-y_i)**2 + (Z-z_i)**2 = (v(t_i - T)**2
其中(x_i, y_i, z_i)
为ith
节点坐标,T
为发射时间
我们现在有 4
个方程 4
个未知数。四个节点明显不够。我们 可以 尝试直接解决这个系统,但是考虑到问题的高度非线性性质,这似乎几乎是不可能的(事实上,我已经尝试 很多 直接技术……失败了)。 相反,我们通过考虑所有 i/j
种可能性,从等式 j
中减去等式 i
,将其简化为一个线性问题。我们得到 (n(n-1))/2
=6
形式的方程:
2*(x_j - x_i)*X + 2*(y_j - y_i)*Y + 2*(z_j - z_i)*Z + 2 * v**2 * (t_i - t_j) = v**2 ( t_i**2 - t_j**2) + (x_j**2 + y_j**2 + z_j**2) - (x_i**2 + y_i**2 + z_i**2)
看起来像 Xv_1 + Y_v2 + Z_v3 + T_v4 = b
。我们现在尝试应用标准线性最小二乘法,其中解是 A^T Ax = A^T b
中的矩阵向量 x
。不幸的是,如果您尝试将其输入到任何标准线性最小二乘算法中,它就会阻塞。那么,我们现在该怎么办?
...
信号到达节点 i
的时间由(当然)给出:
sqrt( (X-x_i)**2 + (Y-y_i)**2 + (Z-z_i)**2 ) / v
这个等式意味着到达时间 T
是 0
。如果我们有 T = 0
,我们可以删除矩阵 A
中的 T
列,问题就会大大简化。事实上,NumPy's
linalg.lstsq()
给出了令人惊讶的准确结果。
...
所以,我所做的是通过从每个方程中减去最早的时间来规范化输入时间。然后,我所要做的就是确定每次可以添加的 dt
,以便最小化线性最小二乘法找到的点的平方和误差的残差。
我将某些 dt
的误差定义为通过将 input times
+ dt
馈送到最小二乘算法预测的点的到达时间的平方差,减去输入时间(标准化),对所有 4
个节点求和。
for node, time in nodes, times:
error += ( (sqrt( (X-x_i)**2 + (Y-y_i)**2 + (Z-z_i)**2 ) / v) - time) ** 2
我的问题:
我能够通过使用蛮力来稍微令人满意地做到这一点。我从 dt = 0
开始,然后逐步移动到某个最大迭代次数 或 ,直到达到某个最小 RSS 错误,这就是 dt
我添加到归一化时间以获得解决方案。生成的解决方案非常准确和精确,但速度很慢。
实际上,我希望能够在 真实 时间内解决这个问题,因此更快 的解决方案将是需要。我开始时假设误差函数(即上面定义的 dt
与 error
)将是高度非线性的——顺便说一句,这对我来说很有意义。
因为我没有实际的数学函数,所以我可以自动排除需要微分的方法(例如 Newton-Raphson
)。误差函数将始终为正,因此我可以排除 bisection
等。相反,我尝试了一个简单的近似搜索。不幸的是,那惨败了。然后我尝试了 Tabu 搜索,然后是遗传算法,还有其他几个。他们都惨遭失败。
所以,我决定做一些调查。事实证明,误差函数与 dt 的关系图看起来有点像平方根,只是根据距信号源节点的距离向右移动:
其中dt
为横轴,纵轴有误差
而且,事后看来,当然有!。我定义了包含平方根的误差函数,所以至少对我来说,这似乎是合理的。
怎么办?
那么,我现在的问题是,如何确定误差函数最小值对应的dt
?
我的第一个(非常粗略的)尝试是在误差图上获取一些点(如上),使用 numpy.polyfit
拟合它,然后将结果提供给 numpy.root
。该根对应于 dt
。不幸的是,这也失败了。我尝试用各种 degrees
和各种点进行拟合,直到点数多得离谱,以至于我还不如使用蛮力。
如何确定这个误差函数的最小值对应的dt
?
由于我们处理的是高速(无线电信号),因此结果 精确 和 准确 很重要,因为次要dt
中的差异可能会影响结果点。
我确信我在这里所做的事情中隐藏了一些非常简单的方法,但是忽略其他所有内容,我如何找到 dt
?
我的要求:
- 速度是最重要的
- 我只能在 运行
的环境中访问纯 Python
和 NumPy
编辑:
这是我的代码。诚然,有点乱。在这里,我使用 polyfit
技术。它将为您“模拟”一个来源,并比较结果:
from numpy import poly1d, linspace, set_printoptions, array, linalg, triu_indices, roots, polyfit
from dataclasses import dataclass
from random import randrange
import math
@dataclass
class Vertexer:
receivers: list
# Defaults
c = 299792
# Receivers:
# [x_1, y_1, z_1]
# [x_2, y_2, z_2]
# [x_3, y_3, z_3]
# Solved:
# [x, y, z]
def error(self, dt, times):
solved = self.linear([time + dt for time in times])
error = 0
for time, receiver in zip(times, self.receivers):
error += ((math.sqrt( (solved[0] - receiver[0])**2 +
(solved[1] - receiver[1])**2 +
(solved[2] - receiver[2])**2 ) / c ) - time)**2
return error
def linear(self, times):
X = array(self.receivers)
t = array(times)
x, y, z = X.T
i, j = triu_indices(len(x), 1)
A = 2 * (X[i] - X[j])
b = self.c**2 * (t[j]**2 - t[i]**2) + (X[i]**2).sum(1) - (X[j]**2).sum(1)
solved, residuals, rank, s = linalg.lstsq(A, b, rcond=None)
return(solved)
def find(self, times):
# Normalize times
times = [time - min(times) for time in times]
# Fit the error function
y = []
x = []
dt = 1E-10
for i in range(50000):
x.append(self.error(dt * i, times))
y.append(dt * i)
p = polyfit(array(x), array(y), 2)
r = roots(p)
return(self.linear([time + r for time in times]))
# SIMPLE CODE FOR SIMULATING A SIGNAL
# Pick nodes to be at random locations
x_1 = randrange(10); y_1 = randrange(10); z_1 = randrange(10)
x_2 = randrange(10); y_2 = randrange(10); z_2 = randrange(10)
x_3 = randrange(10); y_3 = randrange(10); z_3 = randrange(10)
x_4 = randrange(10); y_4 = randrange(10); z_4 = randrange(10)
# Pick source to be at random location
x = randrange(1000); y = randrange(1000); z = randrange(1000)
# 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
t_4 = math.sqrt( (x - x_4)**2 + (y - y_4)**2 + (z - z_4)**2 ) / c
print('Actual:', x, y, z)
myVertexer = Vertexer([[x_1, y_1, z_1],[x_2, y_2, z_2],[x_3, y_3, z_3],[x_4, y_4, z_4]])
solution = myVertexer.find([t_1, t_2, t_3, t_4])
print(solution)
班克罗夫特法似乎适用于这个问题?这是一个纯 NumPy 实现。
# Implementation of the Bancroft method, following
# https://gssc.esa.int/navipedia/index.php/Bancroft_Method
M = np.diag([1, 1, 1, -1])
def lorentz_inner(v, w):
return np.sum(v * (w @ M), axis=-1)
B = np.array(
[
[x_1, y_1, z_1, c * t_1],
[x_2, y_2, z_2, c * t_2],
[x_3, y_3, z_3, c * t_3],
[x_4, y_4, z_4, c * t_4],
]
)
one = np.ones(4)
a = 0.5 * lorentz_inner(B, B)
B_inv_one = np.linalg.solve(B, one)
B_inv_a = np.linalg.solve(B, a)
for Lambda in np.roots(
[
lorentz_inner(B_inv_one, B_inv_one),
2 * (lorentz_inner(B_inv_one, B_inv_a) - 1),
lorentz_inner(B_inv_a, B_inv_a),
]
):
x, y, z, c_t = M @ np.linalg.solve(B, Lambda * one + a)
print("Candidate:", x, y, z, c_t / c)
我的回答可能有错误(明显),因为我在今天下午之前没有听过TDOA这个词。请仔细检查方法是否正确。
我找不到解决您原来的问题的解决方案,即找到与最小误差相对应的 dt。我的回答也偏离了除了 numpy 之外没有第三方库必须使用的要求(我使用 Sympy 并且主要使用来自 here 的代码)。但是,我仍然认为如果所有人都对...感兴趣的话,总有一天有人会发现它很有用……就是找到源发射器的 X、Y、Z。此方法也没有考虑可能存在白噪声或错误或地球曲率和其他并发症的现实情况。
您的初始测试条件如下。
from random import randrange
import math
# SIMPLE CODE FOR SIMULATING A SIGNAL
# Pick nodes to be at random locations
x_1 = randrange(10); y_1 = randrange(10); z_1 = randrange(10)
x_2 = randrange(10); y_2 = randrange(10); z_2 = randrange(10)
x_3 = randrange(10); y_3 = randrange(10); z_3 = randrange(10)
x_4 = randrange(10); y_4 = randrange(10); z_4 = randrange(10)
# Pick source to be at random location
x = randrange(1000); y = randrange(1000); z = randrange(1000)
# 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
t_4 = math.sqrt( (x - x_4)**2 + (y - y_4)**2 + (z - z_4)**2 ) / c
print('Actual:', x, y, z)
我的解决方法如下。
import sympy as sym
X,Y,Z = sym.symbols('X,Y,Z', real=True)
f = sym.Eq((x_1 - X)**2 +(y_1 - Y)**2 + (z_1 - Z)**2 , (c*t_1)**2)
g = sym.Eq((x_2 - X)**2 +(y_2 - Y)**2 + (z_2 - Z)**2 , (c*t_2)**2)
h = sym.Eq((x_3 - X)**2 +(y_3 - Y)**2 + (z_3 - Z)**2 , (c*t_3)**2)
i = sym.Eq((x_4 - X)**2 +(y_4 - Y)**2 + (z_4 - Z)**2 , (c*t_4)**2)
print("Solved coordinates are ", sym.solve([f,g,h,i],X,Y,Z))
根据您给出的初始条件打印语句。
Actual: 111 553 110
几乎立即得出的解决方案是
Solved coordinates are [(111.000000000000, 553.000000000000, 110.000000000000)]
如有不妥,再次抱歉。
背景
我一直在努力解决(众所周知的痛苦)到达时间差 (TDoA) 多点定位问题,在 3-dimensions
中并使用 4
节点。如果你不熟悉这个问题,它是确定一些信号源的坐标(X,Y,Z)
,给定n
个节点的坐标,信号到达每个节点的时间,以及速度信号 v
.
我的解决方法如下:
对于每个节点,我们写(X-x_i)**2 + (Y-y_i)**2 + (Z-z_i)**2 = (v(t_i - T)**2
其中(x_i, y_i, z_i)
为ith
节点坐标,T
为发射时间
我们现在有 4
个方程 4
个未知数。四个节点明显不够。我们 可以 尝试直接解决这个系统,但是考虑到问题的高度非线性性质,这似乎几乎是不可能的(事实上,我已经尝试 很多 直接技术……失败了)。 相反,我们通过考虑所有 i/j
种可能性,从等式 j
中减去等式 i
,将其简化为一个线性问题。我们得到 (n(n-1))/2
=6
形式的方程:
2*(x_j - x_i)*X + 2*(y_j - y_i)*Y + 2*(z_j - z_i)*Z + 2 * v**2 * (t_i - t_j) = v**2 ( t_i**2 - t_j**2) + (x_j**2 + y_j**2 + z_j**2) - (x_i**2 + y_i**2 + z_i**2)
看起来像 Xv_1 + Y_v2 + Z_v3 + T_v4 = b
。我们现在尝试应用标准线性最小二乘法,其中解是 A^T Ax = A^T b
中的矩阵向量 x
。不幸的是,如果您尝试将其输入到任何标准线性最小二乘算法中,它就会阻塞。那么,我们现在该怎么办?
...
信号到达节点 i
的时间由(当然)给出:
sqrt( (X-x_i)**2 + (Y-y_i)**2 + (Z-z_i)**2 ) / v
这个等式意味着到达时间 T
是 0
。如果我们有 T = 0
,我们可以删除矩阵 A
中的 T
列,问题就会大大简化。事实上,NumPy's
linalg.lstsq()
给出了令人惊讶的准确结果。
...
所以,我所做的是通过从每个方程中减去最早的时间来规范化输入时间。然后,我所要做的就是确定每次可以添加的 dt
,以便最小化线性最小二乘法找到的点的平方和误差的残差。
我将某些 dt
的误差定义为通过将 input times
+ dt
馈送到最小二乘算法预测的点的到达时间的平方差,减去输入时间(标准化),对所有 4
个节点求和。
for node, time in nodes, times:
error += ( (sqrt( (X-x_i)**2 + (Y-y_i)**2 + (Z-z_i)**2 ) / v) - time) ** 2
我的问题:
我能够通过使用蛮力来稍微令人满意地做到这一点。我从 dt = 0
开始,然后逐步移动到某个最大迭代次数 或 ,直到达到某个最小 RSS 错误,这就是 dt
我添加到归一化时间以获得解决方案。生成的解决方案非常准确和精确,但速度很慢。
实际上,我希望能够在 真实 时间内解决这个问题,因此更快 的解决方案将是需要。我开始时假设误差函数(即上面定义的 dt
与 error
)将是高度非线性的——顺便说一句,这对我来说很有意义。
因为我没有实际的数学函数,所以我可以自动排除需要微分的方法(例如 Newton-Raphson
)。误差函数将始终为正,因此我可以排除 bisection
等。相反,我尝试了一个简单的近似搜索。不幸的是,那惨败了。然后我尝试了 Tabu 搜索,然后是遗传算法,还有其他几个。他们都惨遭失败。
所以,我决定做一些调查。事实证明,误差函数与 dt 的关系图看起来有点像平方根,只是根据距信号源节点的距离向右移动:
其中dt
为横轴,纵轴有误差
而且,事后看来,当然有!。我定义了包含平方根的误差函数,所以至少对我来说,这似乎是合理的。
怎么办?
那么,我现在的问题是,如何确定误差函数最小值对应的dt
?
我的第一个(非常粗略的)尝试是在误差图上获取一些点(如上),使用 numpy.polyfit
拟合它,然后将结果提供给 numpy.root
。该根对应于 dt
。不幸的是,这也失败了。我尝试用各种 degrees
和各种点进行拟合,直到点数多得离谱,以至于我还不如使用蛮力。
如何确定这个误差函数的最小值对应的dt
?
由于我们处理的是高速(无线电信号),因此结果 精确 和 准确 很重要,因为次要dt
中的差异可能会影响结果点。
我确信我在这里所做的事情中隐藏了一些非常简单的方法,但是忽略其他所有内容,我如何找到 dt
?
我的要求:
- 速度是最重要的
- 我只能在 运行 的环境中访问纯
Python
和 NumPy
编辑:
这是我的代码。诚然,有点乱。在这里,我使用 polyfit
技术。它将为您“模拟”一个来源,并比较结果:
from numpy import poly1d, linspace, set_printoptions, array, linalg, triu_indices, roots, polyfit
from dataclasses import dataclass
from random import randrange
import math
@dataclass
class Vertexer:
receivers: list
# Defaults
c = 299792
# Receivers:
# [x_1, y_1, z_1]
# [x_2, y_2, z_2]
# [x_3, y_3, z_3]
# Solved:
# [x, y, z]
def error(self, dt, times):
solved = self.linear([time + dt for time in times])
error = 0
for time, receiver in zip(times, self.receivers):
error += ((math.sqrt( (solved[0] - receiver[0])**2 +
(solved[1] - receiver[1])**2 +
(solved[2] - receiver[2])**2 ) / c ) - time)**2
return error
def linear(self, times):
X = array(self.receivers)
t = array(times)
x, y, z = X.T
i, j = triu_indices(len(x), 1)
A = 2 * (X[i] - X[j])
b = self.c**2 * (t[j]**2 - t[i]**2) + (X[i]**2).sum(1) - (X[j]**2).sum(1)
solved, residuals, rank, s = linalg.lstsq(A, b, rcond=None)
return(solved)
def find(self, times):
# Normalize times
times = [time - min(times) for time in times]
# Fit the error function
y = []
x = []
dt = 1E-10
for i in range(50000):
x.append(self.error(dt * i, times))
y.append(dt * i)
p = polyfit(array(x), array(y), 2)
r = roots(p)
return(self.linear([time + r for time in times]))
# SIMPLE CODE FOR SIMULATING A SIGNAL
# Pick nodes to be at random locations
x_1 = randrange(10); y_1 = randrange(10); z_1 = randrange(10)
x_2 = randrange(10); y_2 = randrange(10); z_2 = randrange(10)
x_3 = randrange(10); y_3 = randrange(10); z_3 = randrange(10)
x_4 = randrange(10); y_4 = randrange(10); z_4 = randrange(10)
# Pick source to be at random location
x = randrange(1000); y = randrange(1000); z = randrange(1000)
# 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
t_4 = math.sqrt( (x - x_4)**2 + (y - y_4)**2 + (z - z_4)**2 ) / c
print('Actual:', x, y, z)
myVertexer = Vertexer([[x_1, y_1, z_1],[x_2, y_2, z_2],[x_3, y_3, z_3],[x_4, y_4, z_4]])
solution = myVertexer.find([t_1, t_2, t_3, t_4])
print(solution)
班克罗夫特法似乎适用于这个问题?这是一个纯 NumPy 实现。
# Implementation of the Bancroft method, following
# https://gssc.esa.int/navipedia/index.php/Bancroft_Method
M = np.diag([1, 1, 1, -1])
def lorentz_inner(v, w):
return np.sum(v * (w @ M), axis=-1)
B = np.array(
[
[x_1, y_1, z_1, c * t_1],
[x_2, y_2, z_2, c * t_2],
[x_3, y_3, z_3, c * t_3],
[x_4, y_4, z_4, c * t_4],
]
)
one = np.ones(4)
a = 0.5 * lorentz_inner(B, B)
B_inv_one = np.linalg.solve(B, one)
B_inv_a = np.linalg.solve(B, a)
for Lambda in np.roots(
[
lorentz_inner(B_inv_one, B_inv_one),
2 * (lorentz_inner(B_inv_one, B_inv_a) - 1),
lorentz_inner(B_inv_a, B_inv_a),
]
):
x, y, z, c_t = M @ np.linalg.solve(B, Lambda * one + a)
print("Candidate:", x, y, z, c_t / c)
我的回答可能有错误(明显),因为我在今天下午之前没有听过TDOA这个词。请仔细检查方法是否正确。
我找不到解决您原来的问题的解决方案,即找到与最小误差相对应的 dt。我的回答也偏离了除了 numpy 之外没有第三方库必须使用的要求(我使用 Sympy 并且主要使用来自 here 的代码)。但是,我仍然认为如果所有人都对...感兴趣的话,总有一天有人会发现它很有用……就是找到源发射器的 X、Y、Z。此方法也没有考虑可能存在白噪声或错误或地球曲率和其他并发症的现实情况。
您的初始测试条件如下。
from random import randrange
import math
# SIMPLE CODE FOR SIMULATING A SIGNAL
# Pick nodes to be at random locations
x_1 = randrange(10); y_1 = randrange(10); z_1 = randrange(10)
x_2 = randrange(10); y_2 = randrange(10); z_2 = randrange(10)
x_3 = randrange(10); y_3 = randrange(10); z_3 = randrange(10)
x_4 = randrange(10); y_4 = randrange(10); z_4 = randrange(10)
# Pick source to be at random location
x = randrange(1000); y = randrange(1000); z = randrange(1000)
# 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
t_4 = math.sqrt( (x - x_4)**2 + (y - y_4)**2 + (z - z_4)**2 ) / c
print('Actual:', x, y, z)
我的解决方法如下。
import sympy as sym
X,Y,Z = sym.symbols('X,Y,Z', real=True)
f = sym.Eq((x_1 - X)**2 +(y_1 - Y)**2 + (z_1 - Z)**2 , (c*t_1)**2)
g = sym.Eq((x_2 - X)**2 +(y_2 - Y)**2 + (z_2 - Z)**2 , (c*t_2)**2)
h = sym.Eq((x_3 - X)**2 +(y_3 - Y)**2 + (z_3 - Z)**2 , (c*t_3)**2)
i = sym.Eq((x_4 - X)**2 +(y_4 - Y)**2 + (z_4 - Z)**2 , (c*t_4)**2)
print("Solved coordinates are ", sym.solve([f,g,h,i],X,Y,Z))
根据您给出的初始条件打印语句。
Actual: 111 553 110
几乎立即得出的解决方案是
Solved coordinates are [(111.000000000000, 553.000000000000, 110.000000000000)]
如有不妥,再次抱歉。