Более эффективный и надежный алгоритм минимизации цели со многими локальными минимумами. ⇐ Python

[Anonymous]

Общая информация
Общая информация
Я разработал скрипт Python, который вычисляет время, необходимое лучу света для распространения между двумя точками в сложной среде. Теперь я хочу расширить этот сценарий для решения проблемы Разницы во времени прибытия (TDoA) в этой среде. Постановка задачи следующая:
Зная известное расположение N приемников (например, камер, датчиков, антенн) и время прихода светового сигнала на каждый приемник, я стремлюсь определить координаты излучения сигнала. Поскольку время излучения сигнала неизвестно, мне известны только попарные различия во времени поступления (например, разница во времени поступления сигнала между приемником A и приемником B).

Что я уже сделал
Я реализовал сценарий для решения этой проблемы, как показано ниже:

Код: Выделить всё

from dataclasses import dataclass
import numpy as np
from scipy.optimize import differential_evolution
from ray_tracer import RayTracer
from random import randrange

@dataclass
class Optimizer:
ice_model: np.ndarray # Parameters used by my RayTracer class to model complex medium
antenna_coordinates: np.ndarray

def __post_init__(self):
self.ray_tracer = RayTracer(self.ice_model)

def _objective(self, coords: np.ndarray, measured_times: np.ndarray) -> float:
"""Objective function to minimize the error between measured and predicted transit times."""
x, y, z = coords
predicted_times = np.array([
self.ray_tracer.transit_time([x, y, z], [nx, ny, nz])
for nx, ny, nz in self.antenna_coordinates
])

# Compute time difference matrices
measured_time_diff = np.abs(measured_times[:, None] - measured_times)
predicted_time_diff = np.abs(predicted_times[:, None] - predicted_times)

# Calculate error matrix and fill the diagonal with zeros
error_matrix = np.abs(measured_time_diff - predicted_time_diff)
np.fill_diagonal(error_matrix, 0)

# Return the maximum error
return np.nanmax(error_matrix)

def solve(self, measured_times: np.ndarray) -> np.ndarray:
"""Solve the optimization problem to find coordinates minimizing the objective function."""
bounds = [(-1000, 1000), (-1000, 1000), (-1000, 0)]

# Use differential evolution to minimize the objective function
result = differential_evolution(
self._objective,
bounds=bounds,
args=(measured_times,),
strategy='best1bin',  # Strategy for optimization
maxiter=5000,  # Maximum number of iterations
tol=1e-7,  # Tolerance for convergence
popsize=30,  # Population size
mutation=(0.5, 1),  # Mutation factor
recombination=0.7,  # Recombination constant
disp=True  # Display convergence messages
)

if result.success:
return result.x  # Return the optimized coordinates
else:
raise ValueError("Optimization failed: " + result.message)

# Example usage
if __name__ == "__main__":
# Define the ice model
ice_model = np.array([1.78, 0.454, 0.0132])

# Initialize RayTracer with the ice model
ray_tracer = RayTracer(ice_model)

# Generate random antenna coordinates
antenna_coordinates = np.array([[randrange(-100, 0) for _ in range(3)] for _ in range(10)])

# Generate random true coordinates and measured transit times
x, y, z = [randrange(-1000, 0) for _ in range(3)]
measured_times = np.array([
np.random.normal(ray_tracer.transit_time([x, y, z], [nx, ny, nz]) + 10, 1e-12)
for nx, ny, nz in antenna_coordinates
])

# Initialize the optimizer
optimizer = Optimizer(ice_model, antenna_coordinates)

# Attempt to find optimal coordinates
try:
optimal_coordinates = optimizer.solve(measured_times)
print("Optimal Coordinates (x, y, z):", optimal_coordinates)
print("True Coordinates (x, y, z):", x, y, z)
except ValueError as e:
print(e)

Пример вывода:

Код: Выделить всё

Optimal Coordinates (x, y, z): [-384.46240934  730.54140202 -612.43716286]
True Coordinates (x, y, z): -704 -408 -611

• Функция _objective вычисляет время прибытия сигнала в каждое местоположение приемника для заданного предположения о координатах излучения, используя класс RayTracer. Я строю матрицы, содержащие попарные разницы во времени прибытия, и возвращаю наибольшую абсолютную ошибку.
  Метод решения использует оптимизатор Different_evolution SciPy, который я выбрал, потому что он был единственным, который подходил к решению эта проблема.
• Основной блок генерирует пример данных, выбирая случайные координаты приемника и источника, используя RayTracer для вычисления времени прибытия и передавая это время. (с гауссовым шумом и смещением 10, чтобы доказать, что время передачи сигнала не имеет значения для алгоритма) решателю.

Проблема и попытки решения
Мне нужно решение, эффективное с точки зрения времени выполнения, поскольку оно (в идеале) в конечном итоге окажется почти приложение реального времени. К сожалению, хотя я пробовал применить почти все оптимизаторы SciPy (тщательно учитывая каждую из их сильных и слабых сторон и параметров), а также некоторые свои собственные конструкции, мне не удалось точно восстановить расположение источника. Это не улучшается ни при использовании дополнительных приемников, ни при уменьшении гауссовского шума.
Целевая функция сильно нелинейна, и я пробовал применять различные методы сглаживания (сглаживание по Гауссу, потери Хубера, фильтры Савицкого–Голея) безрезультатно. Никакая другая разумная целевая функция не страдает от этой проблемы. Рассмотрим следующий график цели в z, созданный на основе данных примера (и с фиксированными истинными значениями x и y):
< img alt="введите описание изображения здесь" src="https://i.sstatic.net/lQw8XyD9.png" />
Код, опубликованный выше, использует дифференциальную эволюцию решатель, довольно надежно восстанавливает координату z (попадая в пределах единицы или около того примерно в 2/3 случаев и в противном случае не приближаясь к результату), но надежно не может точно восстановить координату x< Координаты /code> и y. Мне приходит в голову, что, возможно, было бы полезно написать какой-нибудь алгоритм, который сначала решает координату z (впоследствии сводя задачу к одному в двух измерениях), однако мне не удалось его написать (я пробовал например, выборка нескольких точек из множества срезов z и выбор среза с наименьшей ошибкой по различным показателям, однако это было медленно и ненадежно).

Вопрос: запрос альтернатив
Какие альтернативные алгоритмы я мог бы использовать, которые были бы быстрее, точнее и надежнее? Я бы предпочел алгоритм с минимальными настраиваемыми параметрами, так как хочу, чтобы его можно было легко применять к различным значениям Ice_model, границам и т. д.

< h3>Изменить:
Я изменил скрипт ray_tracer.py, чтобы он был более производительным, а также немного поигрался с решателем Differential_evolution. Мне удалось убедить его надежно и точно восстановить координату z, а зачастую и координату x. По какой-то причине он редко точно восстанавливает координату y.

Код: Выделить всё

from dataclasses import dataclass
import numpy as np
from scipy.optimize import differential_evolution
from ray_tracer import RayTracer
from scipy.ndimage import gaussian_filter1d
from random import randrange

@dataclass
class Optimizer:
ice_model: np.ndarray  # Parameters used by my RayTracer class to model complex medium
antenna_coordinates: np.ndarray

def __post_init__(self):
self.best_solution = None
self.best_objective_value = float('inf')
self.ray_tracer = RayTracer(self.ice_model)

def _gaussian_kernel(size: int, sigma: float) -> np.ndarray:
"""Create a Gaussian kernel."""
kernel = np.fromfunction(
lambda x: (1 / (sigma * np.sqrt(2 * np.pi))) *
np.exp(-0.5 * ((x - (size - 1) / 2) / sigma) ** 2),
(size,)
)
return kernel / np.sum(kernel)  # Normalize the kernel

def _objective(self, coords: np.ndarray, measured_times: np.ndarray) -> float:
"""Objective function to minimize the error between measured and predicted transit times."""
x, y, z = coords
antenna_coords = self.antenna_coordinates

# Use in-place operations for array manipulation
num_antennas = antenna_coords.shape[0]

predicted_times = self.ray_tracer.transit_time(
np.full_like(antenna_coords, coords), antenna_coords
)

# Precompute and reuse arrays for time differences
measured_diff = np.abs(measured_times[:, np.newaxis] - measured_times)
predicted_diff = np.abs(predicted_times[:, np.newaxis] - predicted_times)

# Create a mask for the upper triangular part (excluding the diagonal)
mask = np.triu_indices(measured_diff.shape[0], k=1)

errors = measured_diff[mask] - predicted_diff[mask]
sigma = 2.0
smoothed_errors = gaussian_filter1d(errors, sigma=sigma)
total_error = np.nansum(np.abs(smoothed_errors))

# Update the best solution found
if total_error < self.best_objective_value:
self.best_objective_value = total_error
self.best_solution = coords

return np.log10(total_error + 1)

def solve(self, measured_times: np.ndarray) ->  np.ndarray:
"""Solve the optimization problem to find coordinates minimizing the objective function."""
bounds = [(-1000, 1000), (-1000, 1000), (-1000, 0)]

# Use differential evolution to minimize the objective function
result = differential_evolution(
self._objective,
bounds=bounds,
args=(measured_times,),
strategy='rand1bin',  # Strategy for optimization
maxiter=500,  # Maximum number of iterations
tol=1e-6,  # Tolerance for convergence
popsize=5,  # Population size
mutation=(1.59, 1.99),  # Mutation factor
recombination=0.7,  # Recombination constant
workers=-1,
disp=True  # Display convergence messages
)

# Return the best solution found, even if maxiter was reached
if result.success or self.best_solution is not None:
return self.best_solution, self.best_objective_value
else:
raise ValueError("Optimization failed without finding a valid solution.")

# Example usage
if __name__ == "__main__":
# Define the ice model
ice_model = np.array([1.78, 0.454, 0.0132])

# Initialize RayTracer with the ice model
ray_tracer = RayTracer(ice_model)

# Generate random antenna coordinates
antenna_coordinates = np.array([[randrange(-100, 0) for _ in range(3)] for _ in range(10)])

# Generate random true coordinates and measured transit times
x, y, z = [randrange(-1000, 0) for _ in range(3)]
measured_times = np.array([
np.random.normal(ray_tracer.transit_time(np.array([[x, y, z]]), np.array([[nx, ny, nz]]))[0] + 10, 0)
for nx, ny, nz in antenna_coordinates
])

# Initialize the optimizer
optimizer = Optimizer(ice_model, antenna_coordinates)

# Attempt to find optimal coordinates
try:
optimal_coordinates = optimizer.solve(measured_times)
print("Optimal Coordinates (x, y, z):", optimal_coordinates)
print("True Coordinates (x, y, z):", x, y, z)
except ValueError as e:
print(e)

Я внес следующие изменения:

• Изменил стратегию с best1bin на rand1bin, чтобы способствовать более рандомизированному поиску (борьба с локальными минимумами).
• Уменьшен размер популяции с 5 до 30. Как это ни странно, меньшие размеры популяции, похоже, дают лучшие результаты.
• Введено распараллеливание, установив worker=-1.
• Изменена функция ошибок к слегка «более гладкой» логарифмической форме.
• Введено небольшое сглаживание по Гауссу.
• Отслеживайте найденное лучшее решение, чтобы даже если решателю «не удается» оптимизировать согласно своему критерию, возвращается что-то (часто это что-то является довольно хорошим, но не отличным результатом).

Подробнее здесь: https://stackoverflow.com/questions/790 ... -many-loca