CC BY
57
ВЕСТН. МОСК. УН-ТА. СЕР. 15. ВЫЧИСЛ. МАТЕМ. И КИБЕРН. 2023. № 3. С. 10-22 Lomonosov Computational Mathematics and Cybernetics Journal
УДК 519.2
Р. А. Гиргидов1
ИМИТАЦИОННОЕ МОДЕЛИРОВАНИЕ РОЯ С ИСПОЛЬЗОВАНИЕМ ПРОИЗВЕДЕНИЯ АДАМАРА
Использование большого количества дронов ставит перед нами задачи корректного, прозрачного и высокопроизводительного моделирования их совместной деятельности в рамках роевого поведения. Для этого была разработан и приведен пример перехода от описания системы управления единичного дрона к поведению роя одинаковых дронов на основе алгебры (произведения) Адамара. Было выяснено, что использование подобного подхода дает повышение производительности модели в десятки раз и позволяет производить моделирование роевых систем с очень большим количеством особей.
Ключевые слова: произведение Адамара, многоагентная система.
Б01: 10.55959/М8и/0137-0782-15-2023-47-3-10-22
1. Введение. Использование дронов в современной жизни приобретает все больший масштаб и в целом перешло из разряда исследований в прикладную область, но исследования по автономным дронам и, в частности, по дронам, объединенным в рой, не только активны, но со временем приобретают все большую значимость. Происходит это из-за того, что при использовании большого количества дронов на передний план начинают выходить не столько проблемы телеуправления и компоновки отдельных дронов, сколько проблемы взаимодействия дронов друг с другом, особенно осуществляемые в автономном режиме. Таким образом, взаимодействие дро-нов внутри роя является важным частным случаем автономного взаимодействия роботизированных систем, при этом важно еще и исследование предполагаемых синергетических эффектов, связанных с количеством и качеством взаимодействий между ними [1].
Традиционно, при описании поведения и моделировании роя используют многоагентную имитационную модель. Для более полной классификации приведем здесь некоторые определения.
Имитационной моделью называют логико-математическое описание объекта, которое может быть использовано для экспериментирования на компьютере в целях проектирования, анализа и оценки функционирования объекта.
Как правило, к имитационным моделям относят модели сложных систем, где известно поведение системы или частей системы в определенных условиях, но целиком описание системы невозможно или очень сложно. Объединение таких частных моделей зачастую выявляет неизвестные ранее связи и особенности поведения моделируемого объекта.
Однако необходимо отметить, что термин "имитационное моделирование" практически не используется в англоязычной литературе [2].
В России к имитационным моделям относят: дискретно событийные (например, модели системы массового обслуживания), модели системной динамики (как правило, имеющие явную схему решения и пригодные для визуализации в реальном времени) и агентные (исследующие взаимодействие похожих агентов).
Моделирование роя обычно относят к имитационному моделированию, так как в основе его лежит многоагентное моделирование. Исследования роевого поведения групп дронов заключаются в выявлении глобальных особенностей роя на основе правил и законов, действующих на каждый дрон в отдельности. Однако многоагентный подход имеет ряд серьезных недостатков, связанных, в первую очередь, с эмпиричностью модели и невозможностью теоретически предсказать свойства роя для получения желаемого поведения. Целью данной работы является применение парадигмы моделирования роя не как набора отдельных агентов, а как единого объекта
1 Санкт-Петербургский политехнический университет Петра Великого, ст. преп., e-mail: [email protected]
исследования. В частности, предлагается набор формальных процедур, позволяющих перейти от существующих многоагентных роевых архитектур к имитационным моделям системной динамики. Данные процедуры должны обеспечить отображение поведения одного дрона в рое на модели роя, обеспечивая повышение производительности модели роя в целом.
Хотя многоагентное моделирование дает вполне разумные результаты, такой способ моделирования не достаточно эффективен и не дает возможности лаконично описать работу роя в целом, а также выявить закономерности его поведения. Нами предлагается методика преобразования описания поведения одного дрона в модель роя целиком на базе произведения Адамара.
Мы исходили из того, что для описания системной динамики роя, необходимо использовать матричное описание модели роя. Однако, в данном случае, стандартные операции линейной алгебры над матрицами не позволили нам получить простое и "прозрачное" описание модели роя, а также разработать формальную процедуру перехода от модели одного дрона и его взаимодействия с соседями к общему описанию роя.
Для решения указанной проблемы, предлагается использовать матричные операции, основанные на произведении Адамара, дополнив его операциями над матрицами различных размерностей, вводя некоторые неканонические, но часто применяемые на практике операции над матрицами.
В рамках данной статьи мы не ставим себе целью провести полное исследование алгебры, основанной на произведении Адамара, а воспользуемся только некоторыми прикладными результатами.
2. Основные понятия и обозначения.
2.1. Алгебраические операции. Рассмотрим основные положения произведения Адамара. Пусть матрицы A и B:
bii ... b\m
A - Ln -
a ii
ani
a im
, B - [bii] -
' L ij J mn
ni
bn
(1)
— прямоугольные матрицы одинаковой размерности (т х п). Тогда введем следующие произведения.
Определение 1. Произведением Адамара будем называть поэлементное умножение матриц А о В = ау Ъу.
Определение 2. Единичной матрицей будем считать матрицу Е = (1)у, все элементы которой равны единице.
Расширим данное определение, определяя его для произведения прямоугольной матрицы на матрицу столбец и матрицу строку.
Определение 3. Введем понятие произведения Адамара для столбца. Пусть X = [ж1,..., хп]т — столбец размерности п. Будем говорить, что
(2)
т.е. каждый столбец матрицы А умножается поэлементно на столбец матрицы X. Для данного определения произведения матрицы и столбца необходимым условием является равная размерность столбцов матрицы и столбца, на который умножается матрица.
Определение 4. Произведением Адамара для строки будем называть следующую операцию. По аналогии со столбцом определим умножение на строку. Пусть У = [у 1 ■ ■ ■ ут] — строка размерности т, тогда
A о X - aj xi - aii . . aim о xi - aiixi . • aimxi
ani • . anm xn anixn • . anmxn
Y о A - aij yj - [yi
Vm]
aii
ani
aim
aiiyi
aniVi
aimym
anmym
(3)
a
о
a
Таким образом, при подобных умножениях размерность результирующей матрицы будет совпадать с размерностью исходной матрицы, умноженной на вектор (записанный в виде строки или столбца), и равняться (т х п).
2.2. Функциональные и дифференциальные операции на основе алгебры Адама-
ра. Функция, действующая на матрицу в рамках данной алгебры, может быть представлена в виде / (А)тхп = [/ (ац)]тхп, т.е. функция действует на каждый элемент матрицы. Символом ^ будем обозначать гиперфункцию ^ (А) = [/ц (ац)], т.е. на каждый элемент матрицы ац действует своя функция /ц. Таким образом, обычная функция, действующая сразу на всю матрицу, будет частным случаем гиперфункции, где все /ц (ац) = / (ац), т.е. одинаковы. Также мы будем пользоваться определением дифференциала функции из [3]. На основании этого определения производная по времени будет рассчитываться следующим образом:
~Ж
ч
Дифференциал матрицы будем задавать как (Ш = [с/ц-].
Определение 5. По аналогии с определением функции введем определение обратной матрицы и опишем обратную матрицу как функцию возведения в степень: А-1 = (а-1)ц.
Как видно из формулы, матрица А не должна иметь нулевых элементов, однако это не очень удобно при практическом использовании, поэтому нами предлагается несколько изменить данное определение, расширив его на матрицы с нулевыми элементами:
А-1 =
х-1
0,
ац = 0,
ац = 0.
(4)
Определение 6. Назовем данную операцию расширенным матричным делением, а матрицу, полученную таким образом, "псевдообратной".
Все приведенные операции справедливы для любых матриц размерностью т х п. Умножение на константу будем рассматривать как частный случай функции, действующей на матрицу.
На основании данных выше определений можно сделать следующие выводы: адамарово произведение матриц одинаковой размерности и матрицы, умноженной на столбец или строку соответствующих размерностей, коммутативно и дистрибутивно;
Апт о Хп + Впт о Хп — (Апт + Впт) о Хп;
(5)
А о (Х^) + Хп2>) = А о Х^) + А о Х^.
(6)
3. Специальные матричные операторы. Введем также оператор понижения размерности матриц. Данную операцию часто называют сверткой. Определение 7. Определим свертку по столбцу как
(А)- = £'
г=1
п
Е аг1
г=1
п
агт
г=1
(7)
где горизонтальная линия в индексе указывает, что результат будет в виде строки, т.е. (А)- — матрица-строка, где каждый ее элемент равен сумме элементов, соответствующего столбца матрицы.
Определение 8. Аналогично введем свертку по строкам как
(А)| = Е'
Н]
3=1
Е аЦ 0 = 1
т
апц
Ь=1
(8)
где вертикальная линия в индексе указывает, что результат будет в виде столбца, т.е. (А)| — матрица-столбец, где каждый ее элемент равен сумме элементов соответствующей строки матрицы А.
Данные соотношения выражаются через стандартное матричное умножение [4] следующим образом: (А)- = Е- ■ А, Е- = [1,1,..., 1], (А)| = А ■ Е|, Е| = [1,1,..., 1]т, где Е| — матрица-столбец, Е- — матрица-строка.
Для корректного описания роевого поведения в случаях различных характеристик дронов в рамках одного роя, описания правил переключения системы управления или логических условий взаимодействия дронов необходимо ввести матричный аналог индикаторной функции.
Для этого зададим гиперфункцию
где Ь С Ф — выбранное подмножество произвольного множества Ф. В матричном виде данная функция будет представлена как
Для дальнейшего описания методики перехода от правил взаимодействия двух дронов к описанию роя на основании введенных определений необходимо будет рассмотреть модель одного дрона и его взаимодействия с другими дронами.
4. Описание модели дрона. Далее в работе рассматривается квадро- или октокоптер с максимальной скоростью полета около 12 м/с. Ограничение скорости в рассматриваемой ниже модели осуществляется двумя способами: прямое ограничение (на практике модель никогда не доходила до срабатывания этого ограничения) и косвенное ограничение (искусственно введенная виртуальная сила вязкого трения). Принципиальной особенностью данного вида дрона (и соответственно модели) является его осесимметричность.
Дрон оборудован оптическими датчиками, позволяющими на близком расстоянии (до та2) определять скорость и координаты дрона соседа или соседей. Их координаты в данной области определяются с некоторой погрешностью (до ±5 см), с частотой, достаточной для определения относительной скорости каждого дрона-соседа. Это зона уверенной детекции. На больших расстояниях (до т) мы можем определить только наличие (отсутствие) дронов в секторе. Конструктивно детекция в секторе представляется фотоэлементом или ультразвуковым датчиком относительно узкоспектрального излучения, направленного в одну сторону и определяющего наличие и интенсивность сигнала от другого дрона. Каждый дрон излучает сигнал в соответствующем спектре. Причем, если дронов в области видимости датчика несколько, то сигнал от них "суммируется" в датчике. Таким образом, получается, что излучение от нескольких далеко стоящих дронов датчиком будет восприниматься как излучение одного, но находящегося ближе и, соответственно, более яркого излучателя. В приведенном ниже примере дрон будет иметь 4 датчика, направленных в разные стороны с растровым углом обзора 90°, модель будет иметь 4 непересекающихся круговых сектора.
Предполагается, что дроны имеют автопилот, благодаря которому все дроны роя находятся на одной высоте. Автопилот присутствует, но может только стабилизировать дрон в горизонтальном положении в отсутствие воздействий от различных управляющих сил и зафиксировать высоту полета дрона. Обратим особое внимание на то, что стабилизация дрона не предполагает остановки дрона, так как у автопилота нет информации о начальной скорости и положении дрона в пространстве и всегда присутствует "паразитный" дрейф.
Движение дрона как материальной точки будет происходить на основании закона Ньютона в дифференциальной форме под действием квазифизических управляющих сил. В данном примере мы будем рассматривать плоскую задачу в предположении, что всем автопилотам роя дали команду держать определенную высоту. В этом случае движение дрона определяется уравнением
1 ь (а)
1, а € Ь, 0, а / Ь,
а / Ь,
(9)
(2 Г -)-)-)-> = + исгг + Рге1 + Ргп(1,
где г = (у) — вектор координат дрона в абсолютной системе отсчета, т — масса дрона, Ffг — результирующая сила взаимодействия одного дрона с другим. Данная сила, как и сила трения, относится к внутренним управляющим воздействиям и зависит от взаимного расположения дронов. Пусть и^г — вектор сил, получаемых от системы управления при обработке внешнего управляющего (целеполагающего) сигнала, ответственного за выполнение задачи (миссии), в отличие от внутренних управляющих сигналов, ответственных за взаимную навигацию дронов. Таким образом, все остальные силы и управляющие команды являются результатом срабатывания правил системы управления при изменении состояния окружающей среды или состояния дрона. Пусть Fгnd — случайное воздействие на дрон, оно может быть как искусственно созданной силой с помощью управляющего модуля, так и силой естественного происхождения, например от порывистого ветрового воздействия или "паразитных" сил, возникающих из-за особенностей работы (неточности) датчиков, а также запаздывания сигнала. Введение данной составляющей управляющих сигналов позволяет добиться получения наиболее устойчивой формации роя, толерантной к большинству внешних воздействий. В данном случае мы используем белошумное воздействие.
Начальные условия задаются как координаты и скорость дрона в начальный момент времени:
го
Жо\ ^ = /Ухо' чУо/' 0 \Vy0y
Рассмотрим силы, действующие на дрон, более подробно. Равнодействующая Ffг от управляющих сил, ответственных за гашение колебаний, моделируется по аналогии с силой вязкого трения и в этом случае будет иметь вид
Ffг — ^ctгV^fг,
где "fг мы рассматриваем как равнодействующий вектор скоростей относительно дрона или дронов, находящихся в зоне уверенной детекции, ^г — управляющий "коэффициент трения", который подбирается экспериментально и зависит от требований, предъявляемых системе управления, связанной с взаимной навигацией дронов.
В случае множества дронов, влияющих на текущий дрон, силу взаимодействия между дро-нами можно представить как сумму сил, создаваемых всеми дронами, находящимися в зоне детекции:
N ^
—гв1 ^ ^ —гв^ • г
Здесь N(¿) — число дронов, находящихся в зоне детекции исследуемого объекта в каждый момент времени ¿. Если рассмотреть описание взаимодействия двух дронов, в данном случае влияния одного дрона на рассматриваемый, то N = 1.
Как было указано выше, все силы, воздействующие на дрон, являются квазифизическими, будучи по сути управляющими воздействиями, но эти воздействия в рамках представленной модели копируют физический процесс — взаимодействие атомов в процессе синтеза кристалла (формирования кристаллической решетки). Это позволяет надеяться на выстраивание дронов в рое в определенную структуру (определяемую параметрами указанных сил).
Отметим важную особенность работы дрона, основанную на физических ограничениях устройства: результирующая сила не может быть больше некой максимальной величины —тах, определяемой индивидуально для каждой модели дрона. В этом случае величина силы становится постоянной и равной максимуму, а вот направление воздействия должно совпадать с направлением результирующей всех управляющих сил, действующих на дрон:
-х =
х гвв
F
F
гв
< —
< - тах) > —тах.
По аналогии с максимальной силой в процессе решения дифференциального уравнения вводится максимально возможная скорость:
Vх =
у твв
V«
Ут
V-,
У«
V,
^ Утах) > Утах)
(11)
но в рамках нашей модели воздействие силы вязкого трения таково, что введение ограничения по скорости практически не используется.
Данные особенности модели дрона описаны в работе [5].
Таким образом, в абсолютных координатах уравнение будет выглядеть так:
= Р/гХ + исЬгХ + КеЯХ + КайХ, = + УсЬтУ + ^ге/У + ^пкй'.
(12)
Подводя итог, выделим некоторые существенные ограничения рассматриваемой нами модели.
• Мы рассматриваем дрон как материальную точку, т.е. отсутствует учет моментов и собственного вращения дрона.
• Мы предполагаем, что дрон может двигаться в произвольном направлении с одинаковыми динамическими характеристиками.
• Мы разделяем внешнее управляющее воздействие, необходимое для целеполагающих задач, стоящих перед дроном (роем) и внутреннего взаимодействия, обеспечивающего взаимную навигацию дронов. При этом все управляющие команды (всех видов) выражаются в виде некоторых виртуальных сил, действующих на дрон, а сам дрон двигается под действием результирующей от всех воздействий.
• Воздействие сил взаимодействия определяется внутренней программой дрона и зависит от используемых датчиков, программно-аппаратных методов идентификации соседей и динамических характеристик дрона.
• Взаимодействующие дроны находятся на одной высоте.
Рассмотрим правила взаимодействия двух дронов более подробно. На рис. 1 выделены следующие области.
Рис. 1. Области видимости датчиков дрона
Светло-серым цветом бозначена область, в которой датчики могут определить наличие дрона с точностью до сектора в 1/4 окружности. Разделение на 4 сектора — это минимально возможное количество секторов, которое дает возможность собрать рой. В этой области действуют силы дальнего притяжения в рамках сектора.
Белым цветом выделен один из секторов, в который попал дрон-сосед (D1). При его нахождении внутри сектора сила взаимодействия определяется аналогично силе дальнего притяжения, действие этой силы определяется особенностью приемных датчиков таким образом, что вектор этой силы направлен вдоль пунктирной линии, проходящей через центр окружности этого сектора. Все секторы равнозначны, т.е. у дрона нет выделенного направления, и если в нескольких секторах есть дроны-соседи, то на дрон действует равнодействующая, пропорциональная яркости излучения в данном секторе (так как используется простейший фотоэлемент, то мы будем считать, что излучение от нескольких дронов в рамках одного сектора будет складываться). Таком образом, один, но находящийся близко дрон будет более значим для системы взаимной навигации, чем несколько, но находящихся на удалении.
Средне-серым цветом обозначена область радиуса ra2, в которой датчики точно определяют расположение и скорость дронов относительно друг друга. В этой области действуют датчики визуального контроля, позволяя определять относительные (относительно положения дрона) координаты соседа. В нашем примере в этой области находится дрон D2. На основании полученной информации мы имеем возможность моделировать силы ближнего притяжения, действующие по направлению к центру масс соседа, а также оценить относительную скорость дрона.
Темнее среднесерого обозначена область радиуса rai, в которой дроны могут находиться неограниченно долго и не взаимодействовать друг с другом.
Темно-серым цветом обозначена область радиуса rd, характеристики датчиков в которой аналогичны описанным для зоны, обозначенной серым цветом. В эту область заходить соседним дронам не рекомендуется, так как это чревато столкновением (так называемая зона безопасности). Внутри этой зоны действуют силы отталкивания.
Стрелками обозначены виртуальные управляющие силы взаимодействия, выделенная стрелка — их результирующая.
Для определения значений и характера сил при взаимодействии дронов будем использовать аналог потенциала Леннарда-Джонса, используемого для описания межмолекулярного взаимодействия [6]. Так, согласно описанию потенциала Леннарда-Джонса, области притягивания и отталкивания определены через воздействия соответствующих сил:
Frel = Fd + Fa + F. Их области воздействия показаны на рис. 2.
Сила отталкивания (desire force) дронов друг от друга согласно аналогии будет задаваться следующим образом:
Fd = кЛо ld (f8) ,
где kd — постоянный коэффициент взаимодействия при отталкивании, r — вектор от текущего дрона как точки отсчета до соседа, взаимодействие с которым мы рассматриваем. Оба дрона в данной модели — материальные точки. Приведем индикаторную функцию расстояния id, определяющую наличие силы отталкивания дронов друг от друга, если они сильно приблизились:
id (r2) = (1 < г2 (13)
10, r2 ^ rd.
Рис. 2. Управляющие силы в зависимости от расстояния до соседа:
F — сила, отвечающая за притяжение дрона на дальних расстояниях;
Fa — сила, отвечающая за притяжение дрона на средних расстояниях;
0 — область отсутствия взаимодействия на ближне-средних расстояниях;
Fd — сила отвечающая за отталкивание дронов друг от друга на ближних расстояниях
Другими словами, сработали датчики сближения или, в другом варианте конфигурации дрона, правила предупреждения столкновений.
Сила Fd виртуальна, т.е. она является результатом работы специального программно-аппаратного модуля, отвечающего за взаимодействие дронов. Она предназначена для двух целей:
1) избегание столкновений;
2) определение минимально допустимого сближения дронов внутри роя.
В общем случае, эти две цели могут быть реализованы разными правилами, но в данном случае действует одно правило.
Необходимо отметить, что на баланс притягивания и отталкивания в случае работы Леннарда-Джонса влияет разность степеней зависимости от расстояния, которая должна быть не менее 2, чтобы у дронов не было возможности с "разгона" при притягивании пролететь всю область безопасности и врезаться друг в друга. Сила притяжения дронов (attraction force) друг к другу на "среднем" расстоянии:
Здесь ка — постоянный коэффициент взаимодействия при притягивании, г — вектор координат соседа, где текущий дрон является точкой отсчета, взаимодействие с которым мы рассматриваем,
1а — индикаторная функция, определяющая наличие силы взаимодействия при притяжении дро-нов друг к другу, если они находятся на определенном расстоянии, эта функция аналогична функции используемой при описании силы "отталкивания", только она действует на средних расстояниях. Данная сила предназначена для обеспечения структурной сборки роя.
Сила притяжения дронов (long distance force) друг к другу на "большом расстоянии":
(14)
(15)
F = ktrti (r2) ,
где к, — постоянный коэффициент взаимодействия при дальнем притягивании, г — единичный вектор определенного углового сектора, сектор определяется наличием фотоэлемента, направленного в определенную сторону, чувствительный элемент которого может "поймать" сигнал в угловом секторе размером п/2 (четверть круга) и определить "светимость". В силу физических особенностей данного датчика сигнал от него пропорционален сумме "светимостей" от нескольких дронов. Необходимо отметить также, что размер сектора выбран таким образом, чтобы обеспечивать приемлемую взаимную навигацию при минимальной стоимости. В процессе работы были также опробованы варианты с размером сектора п/4 и п/6, но существенного выигрыша при стабилизации роя не было.
Приведем индикаторную функцию, определяющую наличие силы взаимодействия при притяжении дронов на дальних расстояниях в рамках сектора,
1,(г»)Л1-
10, г2 ^ г» или
0 0 0 0 (16) 0, г2 ^ г» или г2 ^ га2.
Если дроны совсем далеко, вне радиуса г,, то они друг друга не "видят", а если они оказались внутри радиуса, то максимум, что можно определить — это то, что в данном секторе наблюдается некая "светимость", никакой информации ни о расстоянии до соседей, ни об их относительных координатах, дрон определить не может.
Таким образом, при взаимодействии двух дронов в каждый момент времени на дрон может действовать только одна внутренняя сила в зависимости от расстояния между дронами. Если дронов несколько, то силы их воздействия складываются (а часто и взаимно компенсируются).
В то же время данная модель позволяет снять некоторые ограничения, обычно налагаемые на подобные системы.
1. В нашей модели мы исходим из предположения, что нам не доступна координатная навигация, а дрон с разной степенью точности "видит" только своих соседей. В соответствии с ними работает программа внутреннего управления, подающая команды на контроллер соответствующих двигателей, воздействие которых можно свести к воздействию управляющих сил.
2. Данная модель роя не ограничена законами сохранения (энергии и импульса), однако симметрия пространства и симметрия законов управления дает нам возможность говорить о сохранении положения центра масс роя в отсутствии .
5. Описание состояние роя на основе динамического имитационного моделирования. На основании данной модели самого дрона, модели датчиков и модели управляющих сил (правил) мы построим многоагентную систему роя, где каждый агент взаимодействует с другими через модель окружающей среды.
Запишем состояние роя с использованием матричного подхода. Пусть матрицы
(17)
задают соответствующие координаты каждого дрона в рое в статической системе кординат в каждый момент времени. Тогда Ё (¿) = [X (¿); У (£)] будем называть гипервектором координат роя. Таким образом, попарные относительные координаты дронов будут представлять собой квадратную матрицу:
X = Е о X - Xт о Е, - т (18) У = Е о У - Ут о Е,
где Е = (1)пхп — единичная матрица в смысле произведения Адамара. Здесь столбцы X^, У — координаты группы относительно г-го дрона Приведем пример матриц отклонения X, У для
Х1 (¿) У1 (*)"
X (¿) = ,У (*) =
Хп (¿)_ Уп (*)_
случая трех дронов:
О ж1 - Ж2 ж1 - Ж3 "1 1 1 ж1 1 1 1
Ж2 - ж1 О Ж2 - Жз = 1 1 1 О Ж2 - [xi Ж2 жз] 1 1 1
Ж3 - ж1 Жз - Ж2 О _ 1 1 1 .Жз. 1 1 1
О У1 - У2 У! - Уз "1 1 1 У1 1 1 1
У2 - yi О У2 - Уз = 1 1 1 О У2 - [у! У2 Уз] 1 1 1
Уз - yi У3 - У2 О 1 1 1 _Уз_ 1 1 1
(19)
R2 = X2 + Y2. Для случая трех дронов
О r22 r^'
r2i Lr3i
' 32
2
'23
О
О
(X2 - Xi)2
_(Ж3 - Xi)2 (Ж3 - X2)2
(xi - Ж2)2 (Ж1 - Жз)2 О (Ж2 - Жз)2
О
+
О
(У2 - yi)2
,(Уз - yi)2 (Уз - У2)2
(yi - У2)2 (yi - Уз)2 О (У2 - Уз)2 О
(20)
Введя булеву функцию по аналогии с индикаторной функцией, описанной выше, получим
1гЬг2(я2) = Ь ,. 2 (21)
'2 < r2j < r|,
Г2, ^ или r2- ^ Г2
2 lJ i
Здесь Г1,Г2 — граничные радиусы, внутри которых действует одна из сил взаимодействия п < < т2, т^ — это расстояние от дрона 1 до дрона ,]. Таким образом состояние роя в абсолютных координатах будет описываться уравнениями
m
т-
d2X dt2
dt2
= FfrX + FctrX + FrelX + FrndX, = FfrY + FctrY + FrelY + FrndY •
(22)
В силу линейности изменения матрицы абсолютных и относительных координат в каждый момент времени получаем
dt2 dt2 d2Y d2Y
(23)
dt2 dt2 '
Сила взаимодействия Frei такая же, как для одного дрона, но состоит из гипервекторов. По аналогии с заданным выше гипервектором координат роя мы будем говорить о гипервекторе сил, действующих на рой:
Frei = Fd + Fa + Fi. (24)
В проекции на оси координат:
{FrelX = FdX + FaX + FiX, FrelY = FdY + FaY + Fjy ,
= kd^ld (r2)
Fd,X = kdX о R-5 О Id (R2) : Fd,Y = kdY О R-5 О 1d (R2) ,
(25)
(26)
где Я 5 — псевдообратная матрица Я5. При этом квадратные матрицы взаимного влияния дронов Fdx, FdY при отталкивании преобразуются в столбцы
Fd,X = (Fd,X )|, Fd,Y = (Fd,Y )| •
О
2
Аналогично описывается сила притяжения дронов на близком и на дальнем расстоянии. Отметим, что в данном случае все матрицы симметричны и соответственно симметричны их ада-маровы произведения, а свертка по столбцу будет эквивалентна свертке по строке:
Ра = к
а^Да Г
Ра,Х = кД о Ё
РаХ = каУ о Ё
-3
-3
'1а (Ё2) 1а (Ё2)
(28)
где Ё 3 — псевдообратная матрица Ё3. Аналогично, , Р<1,у — квадратные матрицы взаимного влияния, которые при отталкивании дронов преобразуются к виду
Ра,Х = <Ра,Х)
Ра,Х = <Ра,У)
(29)
6. Выводы. Сравнение двух подходов. Как видно, описание движения единичного дрона при использовании произведения Адамара практически эквивалентно описанию движения всего роя, что очень упрощает как реализацию модели, так и эксперименты с воздействием дронов друг на друга.
С другой стороны, использование такого подхода может накладывать, если не ограничения, то, по крайней мере, может осложнить моделирование в следующих случаях:
• если описание законов управления невозможно свести к воздействию разнонаправленных сил на дрон в каждый момент времени;
• если воздействие дронов друг на друга не является одинаковым, т.е. в рое существуют выделенные (главные) дроны. Данный аспект не исследовался и находится вне рамок данной работы, но наш подход все же позволяет моделировать иерархические модели дронов, если заменить индикаторную функцию или единичную матрицу в формулах относительных координат на весовые матрицы;
• однако самым важным ограничением будет введение таких законов взаимодействия между дронами, при которых воздействие нескольких дронов на один, не аддитивно, однако и в этом случае изменение введенной ранее операции свертки на другую операцию, действующую в рамках описанной алгебры даст возможность обойти это ограничение.
При соблюдении указанных нами ограничений мы можем наблюдать полную идентичность результатов, полученных в расчетах с помощью многоагентной и системно-динамической моделей (см. рис. 3, 4).
■ш
9 ® *
® в ®
в
-4
-2
Рис. 3. Динамическая модель
Рис. 4. Многоагентная модель
о
На рис. 3,4 показано состояние роя через 1 мин 10 с после начала моделирования. Шаг по времени составляет 0.05 с. Начальные состояния совпадают.
Таким образом поставленная нами цель, а именно, получение описания движения роя как единого объекта была достигнута. В процессе разработки модели мы получили методику преобразования частной модели дрона в динамическую модель роя, что было продемонстрировано на приведенном примере перехода от модели одиночного дрона к модели на базе алгебры Адамара.
Одним из важных технических достоинств данного подхода является возможность использовать готовые реализации указанных выше операций на базе произведения Адамара и введенных нами расширений, так как они являются стандартом операций машинного обучения в большинстве современных математических пакетов.
Для сравнения производительности при применении данного подхода нами была реализована описанная выше модель роя на базе многоагентного взаимодействия, а также в виде динамической (матричной) имитационной модели, причем последняя была протестирована в однопоточ-ном, многопоточном вариантах, а также варианте с использованием GPU.
Характеристики решателя для всех реализаций были следующие: состояние роя рассчитывалось в течение внутренних 200 с с шагом 0.05 с, т.е. было сделано 4000 итераций. Опишем подробнее используемые при вычислениях ресурсы.
MultiAgent — мультиагентное моделирование проводилось на шестиядерном процессоре в режиме использования 6 активных ядер для вычислений. Обмен сообщениями производился в режиме синхронного распределения памяти.
Single Thread CPU — динамическое имитационное моделирование на базе алгебры Адамара. Проводилось на том же компьютере с 6 ядрами, но при использовании одного ядра и стандартных библиотек матричных вычислений.
Multi Thread CPU 6 cores — та же самая динамическая модель, но вычисления производились на базе многопоточного вычислителя TensorflowCPU. Использовались все ядра.
Multi Thread CPU 26 cores — модель аналогична предыдущей, но с использованием более мощного процессора. Данная проверка была сделана для подтверждения предположения, что производительность возрастает линейно от количества ядер процессора.
GPU Nvidia 1070 — динамическая модель, аналогичная предыдущим вариантам, но при вычислениях использовалась реализация вычислителя на базе TensorflowGPU.
Моделирование роя как многоагентной системы показало результат в 6 раз медленнее, чем динамическое моделирование в однопоточном режиме, что может быть критически важно для работы в режиме реального времени.
Сравнительные затраты по времени для модели роя, описанного в качестве примера применения алгебры Адамара (с одинаковыми характеристиками решателя)
Количество дронов Multy Agent, с Single Thread CPU, с Multi Thread CPU 6 cores, с Multi Thread CPU 26 cores, с GPU Nvidia 1070, с
100 2110 313 40.2 20 38
500 9641 (3 ч) 1489 196 48 217
1000 20543 (6 ч) 3230 (1 ч) 423 107 413
Как видно из таблицы, рост производительности при задействовании разного количества ядер почти линейный, что говорит о корректности реализованных алгоритмов с точки зрения обработки матриц. Другими словами, предложенный алгоритм описания модели позволил избавиться
от поэлементного вычисления состояния каждого дрона и перейти к высокоэффективным матричным операциям.
Также мы провели эксперименты по расчету поведения роя при использовании GPU в тех же условиях. К сожалению, производительность относительно старых веерных процессоров Pascal based GPU игрового исполнения, сравнима, в данном случае, с использованием 6-ядерного процессора. При этом использование данного подхода осложняется проблемами с отладкой программного обеспечения под GPU. Однако при использовании более мощных современных GPU, ориентированных на высокопроизводительные вычисления (Titanium edition или Tesla) мы получим выигрыш по сравнению с использованием процессоров. Также были произведены тесты масштабирования, которые показали линейный рост времени моделирования относительно роста количества объектов моделирования (дронов).
К существенным недостаткам данного подхода мы бы отнесли следующую особенность: разработка программы для агента фактически означает разработку программы для устройства, при использовании предложенного подхода программу для дрона необходимо будет написать отдельно, но это, как сказано выше, компенсируется возрастанием производительности при моделировании роя почти в 50 раз.
СПИСОК ЛИТЕРАТУРЫ
1. Magnus J.R., Neudecke H. Matrix differential calculus with applications to simple, Hadamard and Kronecker product // J. of Math. Psychol. 1975. P. 3-228.
2. Maria A. Introduction to modeling and simulation // Proceedings of the 29th Conference on Winter Simulation. USA, 1997. P. 7-14.
3. Davis C . The norm of the Schur product operation // Numerische Mathematik. 1962. P. 343-344.
4. Фаддеев Д.К., Фаддеева В.Н. Вычислительыне методы линейной алгебры. Л.: Наука, 1975.
5. G i r g i d o v R.A., Khakimov T. Self-organization of a group of drones in the conditions of absence of data on the exact location // Algorithms and Solutions Based on Computer Technology. Lecture Notes in Networks and Systems. Vol. 387. Springer, 2022. P. 223-233.
6. Jones J.E. On the determination of molecular fields. I. From variation of the viscocity of a gas with temperature // Proceedings of The Royal Society A: Mathematical, Physical and Engineetin Sciences. Vol. 106. 1924. P. 441-462.
Поступила в редакцию 19.01.23 Одобрена после рецензирования 15.03.23 Принята к публикации 15.03.23
