虽然没有类似的任务经验,我读了一些东西,也尝试了一些东西。
当不熟悉这个话题时,很难理解它,而且我发现的所有资源都有点混乱。
的问候理论还不清楚对我来说:
- 高于凸优化问题的规定的问题; (局部最小=全球最低就意味着获得强大的求解器!)
- 大约有更通用的问题更多的资源(传感网 本地化),其是非凸并在极其复杂的方法已经被开发
- 是您的三边测量法,能够利用> 3点(三边测量与多边测量;至少这段代码似乎并不像它可以意味着:噪声性能差)
这里是一些示例代码有两个方法:
- 答:凸优化:SOCP -Relaxation
- B:非线性规划优化
- 使用scipy.optimize
- 非常完美,我合成实验实施;在嘈杂的情况下甚至有好结果;尽管我们使用的数值分化(自动差异很难在这里使用)
一些附加备注:
- 您例B确实有一些(非常糟糕)噪音或者我认为其他一些问题,因为我的方法完全关闭;而特别是B方法照耀我的合成数据(至少这是我的印象)
代码:
import numpy as np
import cvxpy as cvx
from scipy.spatial.distance import cdist
from scipy.optimize import minimize
np.random.seed(1)
""" Create noise-free (not anymore!) fake-problem """
real_x = np.random.random(size=2) * 3
M, N = 5, 10
NOISE_DISTS = 0.1
pos = np.array([(i,j) for i in range(M) for j in range(N)]) # ugly -> tile/repeat/stack
real_x_stacked = np.vstack([real_x for i in range(pos.shape[0])])
Y = cdist(pos, real_x[np.newaxis])
Y += np.random.normal(size=Y.shape)*NOISE_DISTS # Let's add some noise!
print('-----')
print('PROBLEM')
print('-------')
print('real x: ', real_x)
print('dist mat: ', np.round(Y,3).T)
""" Helper """
def cost(x, Y, pos):
res = np.linalg.norm(pos - x, ord=2, axis=1) - Y.ravel()
return np.linalg.norm(res, 2)
print('cost with real_x (check vs. noisy): ', cost(real_x, Y, pos))
""" SOLVER SOCP """
def solve_socp_relax(pos, Y):
x = cvx.Variable(2)
y = cvx.Variable(pos.shape[0])
fake_stack = [x for i in range(pos.shape[0])] # hacky
objective = cvx.sum_entries(cvx.norm(y - Y))
x_stacked = cvx.reshape(cvx.vstack(*fake_stack), pos.shape[0], 2) # hacky
constraints = [cvx.norm(pos - x_stacked, 2, axis=1) <= y]
problem = cvx.Problem(cvx.Minimize(objective), constraints)
problem.solve(solver=cvx.ECOS, verbose=False)
return x.value.T
""" SOLVER NLP """
def solve_nlp(pos, Y):
sol = minimize(cost, np.zeros(pos.shape[1]), args=(Y, pos), method='BFGS')
# print(sol)
return sol.x
""" TEST """
print('-----')
print('SOLVE')
print('-----')
socp_relax_sol = solve_socp_relax(pos, Y)
print('SOCP RELAX SOL: ', socp_relax_sol)
nlp_sol = solve_nlp(pos, Y)
print('NLP SOL: ', nlp_sol)
输出:
-----
PROBLEM
-------
real x: [ 1.25106601 2.16097348]
dist mat: [[ 2.444 1.599 1.348 1.276 2.399 3.026 4.07 4.973 6.118 6.746
2.143 1.149 0.412 0.766 1.839 2.762 3.851 4.904 5.734 6.958
2.377 1.432 0.856 1.056 1.973 2.843 3.885 4.95 5.818 6.84
2.711 2.015 1.689 1.939 2.426 3.358 4.385 5.22 6.076 6.97
3.422 3.153 2.759 2.81 3.326 4.162 4.734 5.627 6.484 7.336]]
cost with real_x (check vs. noisy): 0.665125233772
-----
SOLVE
-----
SOCP RELAX SOL: [[ 1.95749275 2.00607253]]
NLP SOL: [ 1.23560791 2.16756168]
编辑:进一步使用非线性最小二乘代替mo可以实现加速(尤其是大规模)重新使用一般的NLP方法!我的结果仍然是相同的(如果问题会凸出的话,我们的预期是一样的)。 NLP/NLS之间的时间可以看起来像9比0.5秒!
这是我推荐的方法!
def solve_nls(pos, Y):
def res(x, Y, pos):
return np.linalg.norm(pos - x, ord=2, axis=1) - Y.ravel()
sol = least_squares(res, np.zeros(pos.shape[1]), args=(Y, pos), method='lm')
# print(sol)
return sol.x
尤其是第二方法(NLP)也将运行更大的情况下(cvxpy的开销伤害;这不是SOCP求解器应该具备远更好的缺点!)。
这里是一些输出M, N = 500, 1000
一些更多的噪音:
-----
PROBLEM
-------
real x: [ 12.51066014 21.6097348 ]
dist mat: [[ 24.706 23.573 23.693 ..., 1090.29 1091.216
1090.817]]
cost with real_x (check vs. noisy): 353.354267797
-----
SOLVE
-----
NLP SOL: [ 12.51082419 21.60911561]
used: 5.9552763315495625 # SECONDS
所以在我的实验它的工作原理,但我不会放弃任何全局收敛保证或重建,担保(还缺少一些理论)。
起初,我虽然关于使用宽松SOCP问题的全局最优值作为NLP求解器中的初始点,但我没有找到任何需要它的例子!
有的只是换有趣的视觉效果使用:
M, N = 20, 30
NOISE_DISTS = 0.2
...
import matplotlib.pyplot as plt
plt.imshow(Y.reshape(M, N), cmap='viridis', interpolation='none')
plt.colorbar()
plt.scatter(nlp_sol[1], nlp_sol[0], color='red', s=20)
plt.xlim((0, N))
plt.ylim((0, M))
plt.show()
![enter image description here](https://i.stack.imgur.com/2KZB0.png)
而且有些超吵的情况下(不错的表现!):
M, N = 50, 100
NOISE_DISTS = 5
-----
PROBLEM
-------
real x: [ 12.51066014 21.6097348 ]
dist mat: [[ 22.329 18.745 27.588 ..., 94.967 80.034 91.206]]
cost with real_x (check vs. noisy): 354.527196716
-----
SOLVE
-----
NLP SOL: [ 12.44158986 21.50164637]
used: 0.01050068340320306
![enter image description here](https://i.stack.imgur.com/xVZe3.png)
A或B的值是否精确或者它们是否有噪音? –
他们可以被视为相当确切 – marcman
我在问,因为 - 正如sascha也提到的 - B似乎是'错误的'。虽然A显然是正确的,但对于你给定的中心点B应该是类似于[[1.74044707 1.10867308 1.19547313 1.90503438] [1.4004128 0.4014424 0.60096256 1.60036121] [1.70092798 1.04554101 1.13717017 1.86899866]] –