Умный в гору не пойдет29.09.2022 00:16

Доброго дня, хабровчане!

Поделюсь с Вами одним численным алгоритмом вычислительной математики для построения оптимального пути движения на сложной поверхности. Под оптимальностью понимается построение такого маршрута, который проходит по карте местности, огибая все локальные экстремумы (возвышенности и впадины), обеспечивая таким образом некоторую степень минимизации перепада высот и длину всего маршрута.

Входная поверхность задана множеством точек (облаком точек) в трехмерном евклидовом пространстве и имитирует собой ландшафт со всеми его неровностями: возвышенности, холмы, горы, овраги, впадины, разломы и т.д. В качестве начального условия выбирается любая точка на поверхности.

Задачу можно рассматривать как продолжение уже имеющегося маршрута. Куда в конце концов упрется наш искомый путь нам неважно, — главное начать, а там как пойдет, чтобы он обогнул неровности. Алгоритм также не гарантирует прохождения маршрута через какие-либо заранее заданные точки из входного множества точек. Гладкость входной поверхности сильно влияет на гладкость и адекватность получаемого решения.

Алгоритм имеет эвристическую природу, поэтому в данной статье не приводятся строгие математические доказательства существования и единственности найденного решения.

Постановка задачи

В декартовой прямоугольной системе координат на плоскости задана равномерная сетка:

$\{ x_i = i \cdot \Delta x, \ \ i = 0, \ldots, N \} \\ \{ y_j = j \cdot \Delta y,\ \ j = 0, \ldots, M \}$

где (x_i, y_j) — узлы сетки; $\Delta x, \Delta y$ — шаги сеток по осии, соответственно. В каждом узле сетки задано значение $z_{i, j} = z (x_i, y_j)$ , представляющее собой высоту карты (ландшафта) местности в рассматриваемой точке. Значения $z_{i, j}$ , образующие исследуемую поверхность, могут быть отрицательными — в таком случае высоту карты следует понимать как глубину относительно некоторой нулевой отметки z=0 .

В качестве начального условия задана точка $\overline{z}=z(x_{\,\overline{N}}, y_{\,\overline{M}})$ , которая может представлять собой крайнюю точку уже известного маршрута. Его необходимо продолжить далее через всю карту к плоскости x=0 , рисунок 1.

Рисунок 1. Пример ландшафта местности — исследуемая поверхность, заданная множеством точек. Синим цветом отмечена часть известного маршрута. Красная точка — начальное условие, точка из которой необходимо продолжить маршрут.

Численный алгоритм

Рассмотрим три ближайшие к $\overline{z}$ соседние точки поверхности — P_1, P_2, P_3 , находящиеся на соседнем слое по оси . Одна из этих точек (P_2) лежит на том же слое по оси , что и $\overline{z}$ , а остальные (P_1, P_3) лежат по диагонали к точке $\overline{z}$ , рисунок 2.

Рисунок 2. Соседние точки по оси Ox для начального условия. Рисунок 2. Соседние точки по оси Ox для начального условия.

Вычислительный алгоритм построения оптимального маршрута является итерационным. На каждой итерации из трех соседних точек P_1, P_2, P_3 необходимо найти такую, сумма квадрата расстояния и квадрата тангенса угла наклона к которой является минимальной.

Для упрощения выкладок введем следующие обозначения:

$d(\overline{z}, P_k)$ — расстояние от точки $\overline{z}$ до точки P_k, (k =1, 2, 3) ;

$\tan(\overline{z}, P_k)$ — тангенс угла, одна сторона которого — это отрезок, соединяющий $\overline{z}$ и P_k , а другая сторона — отрезок, соединяющий $\overline{z}$ и основание перпендикуляра, опущенного из точки P_k на плоскость перпендикулярную оси и проходящую через точку $\overline{z}$ (углы $\alpha_1, \alpha_2, \alpha_3$ , рисунок 3).

Рисунок 3. Схематический чертеж исследуемого множества точек. Начало искомого пути z, который необходимо продолжить и три соседние точки поверхности — P_1, P_2, P_3, из которых нужно будет вы- брать одну, удовлетворяющую условию минимизации. Точки B_k — это основания перпендикуляров, опущенных из точек P_k

Таким образом, в данных обозначениях смысл вышеописанного алгоритма в следующем: используя начальное условие — последнюю точку известного маршрута, $\overline{z}$ , нам необходимо найти такую точку P_k , чтобы выполнялось условие минимизации:

$\min_{k=1, 2, 3}\left(d^2(\overline{z}, P_k) + \tan^2(\overline{z}, P_k) \right) \tag1$

После нахождения данной точки необходимо принять ее в качестве начального условия для следующей итерации. В результате на каждой итерации мы будем получать новые точки искомого оптимального пути.

Смысл слагаемых в (1) следующий: первое слагаемое отвечает за выбор точки с наименьшим расстоянием до предыдущей, и, как следствие, для минимизации суммарной длины полученного маршрута; второе слагаемое в виде квадрата тангенса необходимо для минимизации перепада высот. Таким образом, вышеописанный алгоритм будет прокладывать маршрут, стараясь обходить все локальные экстремумы, что мы и увидим в
экспериментах. Вместо тангенса можно было бы с тем же успехом использовать синус, что никак не повлияет на результат работы алгоритма.

В процессе расчетов может оказаться, что точка $\overline{z}$ будет находиться на крайнем слое по оси ( y_0 или y_M ) — в таком случае необходимо рассматривать не три соседние точки P_k , а две — одна из которых лежит «по прямой», а другая «по диагонали» к точке $\overline{z}$ , при этом алгоритм построения оптимального маршрута никак не изменится.

Сложность алгоритма составляет $\approx O(3n)$ и зависит от «длины» карты, на которой необходимо проложить искомый маршрут. На каждой итерации необходимо проверять по 2–3 точки и находить минимальное значение (1).

Реализация алгоритма

Исходники находятся здесь

Определим свой класс для точки в пространстве:

class Point3D:
    def __init__(self, N, M, z, PositionalParameter):
        self.N = N
        self.M = M
        self.z = z
        self.PositionalParameter = PositionalParameter

Для удобства вместо и я использовал и , соответственно. PositionalParameter будет равен единице для точек P_1, P_3 и нулю для точки P_2 дабы немного упростить вычисления квадратов, упомянутых в формуле (1).

Определим функцию для подсчета квадрата расстояния:

def GetSquareOfDistance(point, testedPoint, dy):
    """Возвращает относительный квадрат расстояния до проверяемой точки"""
    return math.pow(point.z - testedPoint.z, 2.0) + testedPoint.PositionalParameter * dy * dy

point — это точка $\overline{z}$ , testedPoint — точка P_k , — шаг по оси .

На самом деле здесь не совсем квадрат расстояния, а немного упрощенное значение, которое не учитывает шаг $\Delta x$ , т.к. для поиска минимального значения в (1) его можно выкинуть (предоставлю читателю возможность это проверить).

А также функцию для подсчета квадрата тангенса:

def GetTangentSquare(point, testedPoint, dx, dy):
    """Возвращает квадрат тангенса угла между точками"""
    return math.pow(point.z - testedPoint.z, 2.0) / (dx * dx + testedPoint.PositionalParameter * dy * dy)

В отдельном файле были определены функции для формирования поверхностей как функции двух переменных:

import math

# Возвращает сумму квадрата синуса первого аргумента и квадрата косинуса второго аргумента
# x - первый аргумент
# y - второй аргумент
def SinSquarePlusCosSquare(x, y):
    return math.sin(x) * math.sin(x) + math.cos(y) * math.cos(y)

# Функция Гаусса (двумерная гауссиана). 
# Описание параметров (раздел "Многомерные обобщения"): https://ru.wikipedia.org/wiki/%D0%93%D0%B0%D1%83%D1%81%D1%81%D0%BE%D0%B2%D0%B0_%D1%84%D1%83%D0%BD%D0%BA%D1%86%D0%B8%D1%8F
# x, y              - точка плоскости, в которой необходимо вычислить Гауссиан
# A                 - высота колокола
# sigmaX, sigmaY    - размах колокола по оси Ox и Oy, соответственно
# x0                - сдвиг пика по оси Ox
# y0                - сдвиг пика по оси Oy
def Gaussian(x, y, A, sigmaX, sigmaY, x0, y0):
    xx = (x - x0) * (x - x0) / (2.0 * math.pow(sigmaX, 2.0))
    yy = (y - y0) * (y - y0) / (2.0 * math.pow(sigmaY, 2.0))
    return A * math.exp(-(xx + yy))

Гауссиан отлично подходит для имитации отдельно стоящей горы. А с помощью композиции можно создавать горные цепи разной длины и сложности. SinSquarePlusCosSquare создает, как Вы уже догадались, подобие холмистого рельефа.

В конце расчета маршрут записывается в файл как массив трехмерных точек.

Результаты численного решения

На следующих рисунках показаны результаты численных экспериментов нахождения оптимальных маршрутов. Для имитации ландшафта местности использовались различные аналитические поверхности. Все размерности даны в условных единицах, шаг сеток во всех расчетах по обеим осям был взят равным $\Delta x =\Delta y = 0.1$