, (34)
. (35)
Для решения задачи оптимизации (34), (35) будем применять метод простой итерации, выбирая набор,i= 1,…,Nв (34) с предыдущей итерации. Поиск минимума в (34) позволяет найти набор точекна следущей итерации, по которым определяем новый набор,и так далее. Пусть-- номер итерации,, тогда задача (34), (35) может быть переписана в виде:
,
.
Задача оптимизации,решается итеративно до тех пор, пока последовательностине сойдутся для каждого, при этом считается, что наборопределяет начальное расположение точек, а набор,соответствующий перечень емкостей среды обитания многоугольников Вороного.
Методика расчета матрицы расстояний
Для подсчета обобщенного расстояния (32), (33) между произвольной парой точек на поверхности Земли,построим в начале траекторию малого фрагмента полной окружности Земли, проходящей через заданную пару точек. Считаем, что полная окружность лежит в плоскости, проходящей через центр Земли. Другими словами, построим линию:, где-- параметр, описывающим длину фрагмента окружности. Считаем, что,и, где-- длина дуги фрагмента полной окружности Земли. Определим два единичных по длине вектора, которые указывают на пару выбранных точек,, тогда
. (36)
Пусть векторединичной длины указывает на произвольную точку полной окружности, проходящей через выбранную пару точек. Понятно, что векторлежит в плоскости, образованной векторами. Учитывая (36), после несложных преобразований, найдем
, (37)
где-- угол между парой векторов. Согласно (37) очевидно, чтои.
Учитывая (37), а также считая, что, найдем параметрическую запись полной окружности в координатах “широта - долгота”:
, (38)
.
Отметим, что формулаверна, когда. В двух других случаях: 1); 2)к выражению в (38ў) необходимо добавить и вычестьpсоответственно.
При проведении линии между парой точек согласно формулам (38),необходимо также найти ее пересечение с береговой линией для выделения тех частей траектории, которые лежат отдельно на суше и на море. Для изображения траекторий (38),необходимо различать два случай: 1); 2). В первом случае нулевой меридиан поместим в центр карты, т.е. в центре будет располагаться Атлантический океан. Во втором случае после сдвига долготыпостроим карту с центром по линии смены дат, т.е. с центром в Тихом океане.
Рис.17. Набор малых фрагментов полных окружностей, соединяющих пары точек
На рис.17 приведен пример позиционирования десяти точек на поверхности суши. Точки выбраны случайно согласно алгоритму раздела №6, т.е. с учетом плотности емкости среды обитания. Точки пронумерованы и изображены на каждой из двух карт. На рис.17 с учетом сферической геометрии построены также вселиний, соединяющих каждую пару точек. Кроме того, линии размечены в части их прохождения по суше (сплошная линия), и по морю (пунктирная линия). Согласно рис.17 некоторая часть бинарных линий позиционирована на правой карте с центром в Тихом океане.
Согласно формуле (32) для подсчета эффективного расстояния по перемещению одной условной единицы массы груза на поверхности Земли важно знать производную рельефапо маршруту, который опишем некоторой траекторией. Пусть, например, между парой точек,проложен маршрут такой, что, где-- параметр, описывающий пройденный путь от начальной точки. В этом случае очевидно, что.
В качестве маршрута между парой точек на поверхности Земной сферы выберем меньший фрагмент полной окружности, проходящей через пару точек. Подходящая траектория представлена в виде формул (38),. Осталось найти производные. После несложных выкладок получим:
, (39)
.
Опишем большой фрагмент полной окружности, проходящей через некоторую пару точек на поверхности Земли. С учетом описания малого фрагмента полной окружности (37), (38),, параметризацию большого фрагмента целесообразно произвести в два этапа: 1) для значений параметра длины дугииз диапазона, где; 2) для значений параметра длины дугииз диапазона. Уголнайдем из условия того, что долготасовпадает с меридианом смены дат, т.е.. Из решения последнего уравнения найдем, когдаи, когда, при этом
.
Рис.18. Разметка полной окружности на три фрагмента
На рис.18 приведен пример полной окружности, проходящей через пару случайно выбранных точек на поверхности Земли. Окружность поделена на три дуги: 1) малый фрагмент полной окружности помечен маркерами в виде звезд; 2) большой фрагмент полной окружности, параметризованный отрезком дугииз диапазона, где, помечен пентаграммами; 3) большой фрагмент полной окружности, параметризованный отрезком дугииз диапазонаи помечен гексаграммами. Каждая из дуг отмечена соответствующим углом. Стрелка, отделяющая пентаграммы и гексаграммы, обозначает линию смены дат. После объединения второго и третьего фрагментов дуг получим большой фрагмент полной окружности.
Введем матрицу расстояний между всеми парами точек.
Определим подобные матрицы расстояний, найденные для малых и большихфрагментов дуг соответствующих полных окружностей.
Учитывая решение задачи минимизации затрат на трафик, положим, что искомая матрица расстоянийявляется поэлементным минимумом пары матриц расстоянийи, т.е.
.
Проведем вычислительный эксперимент по оценке среднего значения транспортных интегралов в формулах (32), (33). Для подсчета транспортных интегралов нам необходимы частные производныеи, которые вычислим с помощью конечных разностей в форме (17), (18) с разрешением рельефа; а также обыкновенные производныеи, которые найдем согласно (39),. В остальном фрагмент полной окружности, соединяющий пару точек, делился на части, проходящие отдельно по суше и по морю. Для каждой части строилась конечно-разностная сетка с числом точек, равной величине, где-- функция целой части числа, а парыиобозначают начало и конец (в градусах), рассматриваемой части. В дальнейшем в ряде случаев в целях экономии вычислительных ресурсов параметр 0,25 заменялся на 0,125. В итоге для каждой пары точек из общего числавычислим транспортные интегралы и найдем средние значения из каждого набора.
Таблица №7.Усредненные значения транспортных интегралов
В таблице №7 приведен итог вычислительного эксперимента в формате, когда случайно согласно алгоритму раздела №6 выбиралосьточек и вычислялась матрица. Черта сверху во второй строке таблице №7 обозначает операцию усреднения. С учетом оценок из предыдущего разделаследует, что каждое слагаемое в (32), (33) с точки зрения средних значений вносит сравнимый вклад.
В итоге остановимся на выборе значений параметров.
Согласно определению в (32), (33) и таблице №7 интегралы берутся по маршруту, соединяющему пару точек на поверхности Земли. Поскольку все маршруты располагаются на одном из фрагментов полных окружностей, транспортные интегралы можно приближенно считать пропорциональными длине дуги. С этой точки зрения в среднем можно полагать, что матрица расстоянийсовпадает с матрицей расстояний, найденной для малых фрагментов полных окружностей, т.е. считаем в дальнейшем, что.
Алгоритм минимизации транспортных издержек
Вернемся к задаче оптимизации транспортных издержек в формес учетом: 1) особенностей геометрии “суша - море” и 2) вычислительной трудоемкости подсчета функционала транспортных издержеки, в особенности, его частных производных по. Пусть на предыдущем шаге итерационной процедуры по формулам (35ў) подсчитаны объемы емкостей среды обитания, что позволяет определить функционал транспортных издержекв рамках текущей итерации.
Для поиска минимума функциирассмотрим несколько модифицированную схему градиентного спуска. Запишем следующую систему обыкновенных дифференциальных уравнений:
, (40)
где-- некоторый вспомогательный аргумент,,;;,-- так называемая “знаковая” функция, а наборы неотрицательных коэффициентовиопределим ниже.
Заменим обыкновенные производные в (40) на конечные разности, тогда
, (41)
где-- конечные приращения функций, а-- конечное приращение аргумента.
С учетом (41) запишем алгоритм перехода от текущих значений функцийпри значении аргументак новым значениям функцийпри значении аргумента, а именно
,
где. Если после пересчета согласно алгоритмуодна или несколько точек выходят за пределы суши, уменьшаем соответствующие коэффициенты из наборовитак, чтобы точки вернулись в пределы суши. В случае, если береговая линия не мешает движению точек, соответствующие значения коэффициентов полагаются равными единице.
Значение параметрав (41),выбиралось из набора шаговвида:
, (42)
где,-- некоторые константы, а-- натуральные значения. Конкретные значения параметровивыбирались при проведении численных расчетов.
С учетом (41ў), (42) подсчитаем новые положения точекпо формулам:
, (43)
с шагами, а также значения функционала транспортных издержек:
. (44)
Пусть в наборе (44) при просмотре слева направо находится номер, при котором транспортные издержки становятся меньше тех, которые имели место на предыдущем шаге схемы градиентного спуска, т.е.
. (45)
В случае выполнения неравенства (45) при некоторомсчитается, что шаг градиентного спуска завершен и новые положения точек считаются равными:
. (46)
Если неравенство (45) не выполняется при всех значениях, то считается, что данный шаг процедуры градиентного спуска также завершен и новые положения точек выбираются равными:
.
Для замыкания алгоритма минимизации (40) --необходимо уточнить процедуру подсчета номерав наборе (42). Обозначим числос предыдущего шага процедуры--символом, тогда номерв наборе (42) подсчитаем по формуле:и,.
Согласно алгоритму (40) --необходимо иметь частные производные функции транспортных издержекпо координатам точек. Поскольку сделать это аналитически затруднительно, имеет смысл рассмотреть конечные разности. Введем некоторые приращения аргументови запишем следующие выражения:
, (47)
где. При этом, как и выше, будем считать, чтоNточек с координатаминаходятся в пределах суши. Если это не так варьируем приращения. С учетом (47) положим, что
. (48)
Выполняем процедуру градиентного спуска (40) -- (48) до тех пор, пока интегральные издержки на трафик понижаются. Если дальнейшее понижение трафика не удается, считаем согласно формуле, что найдены новые положения Nточек. Пересчитываем по новым значениям Nточек объемы емкостей среды обитаниясоответствующих многоугольников Вороного. Переходим к следующему шагу, т.е. к оптимизационной процедуре (40) -- (48).
Невозможность дальнейшего понижения трафика означает, что шаг градиентного спускастановится минимальным и равным, а доля вариации затрат на трафик за один шаг градиентного спуска периодически меняет знак так, что интегральные затраты на трафик колеблются возле некоторого постоянного значения.
В процедуре градиентного спуска (40) -- (48) с целью понижения вариабельности при конечно-разностном подсчете производных транспортных издержек (47), (48), было проведено сглаживание рельефа с помощью процедуры скользящих средних. Усреднение проводилось по ближайшим соседям в области квадратной формы в координатах “широта - долгота” с числом точек, при этом исходная матрица рельефа имела разрешениеи была представлена в форме матрицы размером.
Для исходной оценки параметра, входящего в формулы, (34),, будем исходить из приближенной оценки доли транспортных затрат в мировом ВВП в диапазоне от 4% до 15% в зависимости от развитости транспортной инфраструктуры. В расчетах будем полагать эту долю, равной приблизительно 10%.
Параметроценим, исходя из формулы, полагая, что затраты на трафик составляют некоторую долюот всей емкости среды обитания, т.е.
. (49)
Считая заданными: константу, набор, матрицуи долюрешим уравнение (50) и найдем коэффициент, который должен обеспечить на финальной стадии процесса минимизации,долю трафика от всей емкости среды обитания в окрестности 10%.
Рис.19,а. Зависимость затрат на транспорт от параметра
Рис.19,б. Многоугольники Вороного для начального распределения cточек
На рис.19,а согласно (49) построен типичный пример графика функциипри. Пентаграммой на графике отмечено значение параметра, которое получено путем решения уравнения (49) при. Матрица расстоянийбыла подсчитана по случайному набору точек, полученному согласно алгоритму раздела №6. На рис.19,б приведены позиции случайного начального распределение точек (обозначены центрами маркеров в виде окружностей o) на поверхности Земли, а также соответствующие многоугольники Вороного.
На рис.20,а,б приведена динамика зависимости транспортных затрат от номера шага градиентного спуска согласно процедуре (40) -- (48) с параметрами:,,для двух наборов точек числом 36 и 72 соответственно.
Рис.20,а. Динамика зависимости транспортных затрат от номера шага градиентного спуска дляточек
Рис.20,б. Динамика зависимости транспортных затрат от номера шага градиентного спускаточек
Для 36 точек в целом было осуществлено 26 шагов процедуры (40) -- (48), затраты на трафик при этом уменьшились на 5,7% по сравнению с исходным значением. Окружностями на рис.20,а,б обозначены значения транспортных затрат после каждого шага процедуры градиентного спуска. Итоговая доля затрат на трафик для 36 точек по отношению ко всей емкости среды обитания составила 10,4%. В течение всего расчета один раз пересчитывался набор емкостей среды обитания. На рис.20,а,б места пересчета обозначены стрелками. На завершающей стадии расчетов на рис.20,а,б виден выход значений трафика на некоторое плато.
На рис.21 приведен итог расчета положенийточек после применения процедуры оптимизации затрат на транспорт (40) -- (48) и одного пересчета набора емкостей среды обитания, ассоциированных с каждой из точек. Первоначально точки выбирались случайно согласно алгоритму (21) -- (24), они помечены центрами маркеров в виде синих окружностей (o) на рис.19,б и на рис.21. После применения процедуры оптимизации,, (40) -- (48) положение точек помечено маркерами в виде пентаграмм. Максимумы смещений составилиипо широте и долготе соответственно.