Tout en n'ayant aucune expérience avec des tâches similaires, j'ai lu quelques trucs et aussi essayé quelque chose.
Lorsque vous ne connaissez pas ce sujet, il est difficile à saisir et toutes les ressources que j'ai trouvées sont un peu chaotiques.
encore peu claire en ce qui concerne la théorie pour moi:
- est le problème comme indiqué ci-dessus un problème d'optimisation convexe; (-minimum local =-minimum global signifierait l'accès aux solveurs puissants!)
- il y a beaucoup plus de ressources sur les problèmes plus génériques (réseau de capteurs Localisation), qui sont non-convexe et où des méthodes extrêmement complexes ont été développés
- est votre approche de trilatération capable d'exploiter> 3 points (trilatération vs multilatération; au moins ce code ne semble pas comme il peut ce qui signifie: mauvaise performance avec le bruit)
Voici quelques exemples de code avec deux approches:
- A: Convex-optimisation: SOCP -relaxation
- B: optimisation Nonlinear programmation
- Mis en œuvre en utilisant scipy.optimize
- Quasiment parfait dans mes expériences synthétiques; même de bons résultats dans un cas bruyant; malgré le fait que nous utilisons la différenciation numérique (disque automatique diff à utiliser ici)
une remarque supplémentaire:
- Votre exemple B a sûrement quelques-uns (très mauvais) bruit ou un autre problème à mon avis, car mes approches sont complètement éteintes; tout particulièrement l'approche B brille pour mes-données synthétiques (au moins c'est mon impression)
code:
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)
Sortie:
-----
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]
Edit: plus l'accélération peut être atteinte (surtout à grande échelle) en utilisant les moindres carrés non linéaires au lieu du mo re approche générale NLP! Mes résultats sont toujours les mêmes (comme prévu si le problème serait convexe). Les délais entre NLP/NLS peuvent ressembler à 9 vs 0.5 secondes!
Ceci est ma méthode recommandée!
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
Surtout exécutera également des instances beaucoup plus la deuxième approche (PNL) (les frais généraux de cvxpy fait mal, ce n'est pas un inconvénient du SOCP résolveur qui devrait évoluer beaucoup mieux!).
Voici quelques sortie pour M, N = 500, 1000 avec un peu plus de bruit:
-----
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
Donc, dans mes expériences, il fonctionne, mais je ne donne aucune garantie globale de convergence ou reconstruction garanties (manque encore une théorie). Dans un premier temps, j'ai pensé à utiliser l'optimum global du problème relaxé-SOCP comme point initial dans le solveur NLP, mais je n'ai trouvé aucun exemple où cela soit nécessaire!
Quelques visuels juste pour-fun en utilisant:
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()
Et certains cas, super bruyant (belle performance!):
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
Les valeurs de A ou B sont-elles des valeurs exactes ou sont-elles bruitées? –
ils peuvent être traités comme assez exact – marcman
Je demande parce que - comme sascha également mentionné - B semble être «faux». Bien que A est évidemment correct, B pour votre point central donné devrait être quelque chose comme [[1,74044707 1,10867308 1,19547313 1,90503438] [1,4004128 0,4014424 0,60096256 1,60036121] [1,70092798 1,04554101 1,13717017 1,86899866]] –