Об оптимизации реализации алгоритма локализации инвариантных множеств динамических систем

Сер. 10. 2010. Вып. 2

ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА

УДК 517.9

С. В. Терентьев

ОБ ОПТИМИЗАЦИИ РЕАЛИЗАЦИИ АЛГОРИТМА ЛОКАЛИЗАЦИИ ИНВАРИАНТНЫХ МНОЖЕСТВ ДИНАМИЧЕСКИХ СИСТЕМ

1. Введение. Компьютерное моделирование динамических систем со сложным поведением является важной частью их исследования. Одна из основных характеристик динамической системы - ее инвариантные множества. Их локализация позволяет выявить динамику системы. Поскольку аналитическое описание таких множеств редко бывает возможным, большое внимание уделяется разработке и реализации компьютерных методов их построения. На практике хорошо зарекомендовали себя методы, основанные на идее аппроксимации фазового пространства конечным набором многомерных кубов (ячеек) [1-3]. Построенное приближение дает представление о поведении системы. С помощью техники последовательного подразбиения ячеек можно уточнить результат. При стремлении диаметров ячеек к нулю и выполнении некоторых ограничений на систему этот процесс сходится и позволяет получить приближение к инвариантному множеству с нужной точностью. Выбор представления ячейки существенно влияет не только на точность результата, но и на скорость работы системы. Мы рассматриваем ячейки разбиения как интервальные вектора в пространстве соответствующей размерности, используем методы интервальной арифметики для построения образа ячейки, а также представление системы с помощью символического образа [1], который является ориентированным графом, созданным по исходной системе.

Реализация алгоритма локализации инвариантных множеств динамических систем (ДС), основанная на применении методов интервальной арифметики, была описана в работе [4] и обрабатывала дискретные ДС размерности 2. Позднее была введена поддержка непрерывных ДС размерности 2 и 3, а также дискретных ДС размерности 3 [5]. В данной статье подробно рассматриваются вопросы оптимизации алгоритма локализации инвариантных множеств ДС. Такая проблема возникла при обработке непрерывных ДС и ДС размерности 3, так как время, требуемое на их обработку, было существенно больше времени обработки дискретных ДС размерности 2. После анализа исходных кодов программы в [5] были предложены техника смешанных вычислений, позволяющая эффективно работать с описанием ДС, и распределенный алгоритм обработки ячеек покрытия.

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

Терентьев Сергей Валерьевич — аспирант кафедры информатики математико-механического факультета Санкт-Петербургского государственного университета. Научный руководитель: кандидат физико-математических наук, доц. Н. Б. Ампилова. Количество опубликованных работ: 7. Научные направления: компьютерное моделирование динамических систем, символическая динамика. Е-mail:waterq<Э.yandex.ru.

© С. В. Терентьев, 2010

с покрытием. В данной работе показано, что при организации хранения ячеек покрытия специальным образом скорость обработки ДС повышается.

2. Локализация инвариантных множеств.

2.1. Постановка задачи. Пусть M - замкнутое и ограниченное подмножество в Rq. Рассмотрим дискретную динамическую систему, порожденную гомеоморфизмом f : MM ^ MM. Пусть C0 = I Mil,  , Mn} - начальное разбиение множества MM.

Символический образ системы f [1] представляет собой ориентированный граф, построенный по системе f и конечному покрытию (разбиению) фазового пространства. Вершины графа соответствуют ячейкам покрытия, а дуга i ^ j строится, если пересечение f (Mi) с Mj не пусто. Более мелкое покрытие позволяет получить более точное приближение к динамике системы. Основные числовые характеристики символического образа (диаметр покрытия, диаметр образа покрытия и нижняя грань) позволяют оценивать точность приближения. Для непрерывных систем символический образ строится по отображению сдвига вдоль траектории.

В работе [4] был рассмотрен алгоритм локализации инвариантных множеств, основанный на представлении системы с помощью символического образа и представлении ячеек покрытия в виде интервальных векторов. При этом мы фактически не строим граф, а работаем с множеством, содержащим ячейки, участвующие в формировании символического образа. При последовательном измельчении разбиения можно получить более точный результат. После осуществления к подразбиений и применений алгоритма результатом вычислений будет последовательность приближений инвариантного множества {Co, C1, Ck}, где

Ci = M,  ,MMjn\Mj1 с Ci-1}, i (={1,.  ,к}.

2.2. Реализация. В настоящее время существуют следующие реализации этого алгоритма локализации инвариантных множеств:

1) базовая (ISLM) [4];

2) модификация ISLM (ISLM+) [5].

Все приложения реализованы на языке C+—Ъ При реализации алгоритма была применена интервальная арифметика библиотеки BOOST [6]. Визуализация инвариантных множеств осуществляется с помощью технологии GNUPLOT [7].

В ISLM была реализована поддержка дискретных систем размерности 2. Реализация ISLM+ является улучшенной модификацией ISLM. В эту программу была внедрена поддержка непрерывных систем размерности 2, дискретных и непрерывных систем размерности 3, а также был улучшен процесс локализации инвариантных множеств. В ISLM+, так же как и в ISLM, в качестве покрытия(хранилища ячеек) применяется структура данных - множество (SET), а операции поиска, вставки и удаления производятся за время, пропорциональное logw, где N - количество элементов в множестве. Моделью данных, соответствующей SET, служит двунаправленный список, все элементы которого различны.

Так как операциями, существенно влияющими на скорость работы ISLM, являются вычисление образа ячейки (Img) и пересечение образа ячейки с покрытием (Intersect), то на их выполнение тратится основная часть общего времени работы программы. Заметим, что процесс создания образа ячейки в ISLM был организован не очень эффективно: при вычислении арифметических выражений необходимо было обойти специально построенное дерево разбора выражений, описывающих искомую динамическую систему. В ISLM+ данную проблему удалось избежать путем внедрения

механизма смешанных вычислений [5], при этом существенно повысилась эффективность операции.

Процесс локализации инвариантных множеств ДС стал более эффективным благодаря не только механизму смешанных вычислений, но и реализации распределенного алгоритма этих множеств с помощью библиотеки Boost.Threads, которая входит в состав Boost - собрание библиотек, расширяющих C++. Текущая реализация такой библиотеки работает на платформах POSIX, Win32 и Macintosh Carbon. При реализации распределенных алгоритмов возникает задача организации доступа к разделяемым ресурсам: необходимо исключить ситуацию, при которой в один момент времени один из потоков меняет значение разделяемых данных, а другой поток пытается его прочитать. В качестве средства синхронизации потоков процесса был использован Мьютех (англ. mutex, от mutual exclusion) - механизм взаимного исключения (рис. 1). Распределенный алгоритм основан на идее распараллеливания операций Img и Intersect.

Рис. 1. Модель потока данных ISLM+

В каждом потоке-обработчике выполняются обе операции. Количество таких потоков варьируется, зависит от аппаратной конфигурации и настраивается опытным путем. Поток-наполнитель всегда присутствует в единственном экземпляре, так как несколько потоков не могут одновременно производить модификации нового покрытия. В данном алгоритме организован синхронный доступ к ячейкам текущего покрытия: недопустима ситуация, при которой одна ячейка будет обработана несколькими потоками, так как производительность программы в таком случае может понизиться. Асинхронный доступ потоков-обработчиков и потока-наполнителя ведет также к возникновению исключительных ситуаций в работе ISLM+.

3. Использование ^-деревьев.Так как операция Intersect сводится к перебору ячеек покрытия, то предлагаемая оптимизация ориентирована на выявление способов организации представления покрытия таким образом, чтобы поиск ячейки становился наиболее эффективным. Был применен подход, предложенный в работе [8]. Он основан на понятии ^-деревьев. В этом случае фазовое пространство разбивается

на последовательность иерархически вложенных многомерных кубов, а процесс поиска ячейки представляет собой выявление и отсечение тех областей покрытия, в которых ячейки быть не может. Модель потока данных в этом случае такая же, как и в ISLM+.

3.1. Основные понятия. Д-дерево - древовидная структура данных (сбалансированное дерево), которая используется для индексации многомерной информации при организации доступа к пространственным данным [8]. Она была предложена в 1984 г. Антонином Гуттманом. Эта структура данных активно применяется в многомерных СУБД для быстрой обработки сложных запросов к БД. Типичным примером использования Д-дерева будет следующая ситуация: найти все банкоматы, магазины, кинотеатры в радиусе 3 м от указанного местоположения (задается широтой и долготой).

Узлы Д-деревьев делятся на листовые и не листовые.

Каждая листовая вершина содержит список индексных записей следующего вида: (I,id), где id - индефикатор, соответствующий многомерному объекту в многомерной базе данных, использующийся для доступа к нему; I = (Io,Ii,..., In-i) - минимальный многомерный куб, содержащий многомерный объект, определяемый индексом id; n -размерность пространства; Ii - замкнутый интервал вида [a, b], a,b G Д.

Каждая не листовая вершина содержит список записей такого вида: (Iu ,child-pointer), где child-pointer - указатель на дочернюю вершину; Iu - минимальный многомерный куб, содержащий многомерные кубы дочерних узлов дочерней вершины.

Для листовых вершин вводится количество индексных записей, для не листовых -число дочерних узлов, а также определяются максимальное (M) и минимальное (m) их значения. При этом то ^ Ц-.

Д-деревья обладают свойствами:

1) количество индексных записей в каждой листовой вершине - от m до M, иначе это корневая вершина;

2) в каждой индексной записи вида (I, id) I - минимальный многомерный куб, содержащий многомерный объект, определяемый индексом id;

3) количество дочерних узлов в каждой не листовой вершине m - M, иначе это корневая вершина;

4) листовая вершина без родительского узла является корневой;

5) в каждой записи вида (Iu,child-pointer) Iu - минимальный многомерный куб, содержащий многомерные кубы дочерних узлов вершины, на которую указывает child-pointer;

6) если у вершины нет дочерних вершин, то она листовая;

7) все листовые вершины находятся на одном уровне.

Максимальная высота дерева, содержащего N индексных записей, равна |logmN| — 1,

1 i N i N

максимальное количество узлов — 1 + — + .

J m m2

На рис. 2 приведена схема расположения узлов Д-дерева в памяти. Д-дерево является обобщенным вариантом Д-дерева с многомерным индексом в качестве ключа.

На рис. 3 показано разбиение многомерного пространства на множество иерархически вложенных многомерных кубов, которые могут пересекаться. Схема связи объектов представляет собой Д-дерево. Поиск индексных записей определяется следующим образом: найти все индексные записи, пересекающиеся с многомерным кубом (запросом). Операция поиска использует ограничивающие кубы для выбора вершин, в которых будет происходить поиск: рассматриваются только вершины, которые пересекаются с запросом. При таком подходе будет осуществляться обход незначительного числа вершин.

Data

Рис. 2. Представление Ä-дерева в памяти

Рис. 3. Логическая структура Я-дерева

Операция добавления комбинируется с операцией поиска: при добавлении индексной записи производится поиск листовой вершины, добавление в которую повлечет

минимальное расширение ограничивающего ее куба. Если после добавления листовая вершина переполнена (то > Ц-), то происходят ее расщепление и перестройка дерева, осуществляемая с помощью многократного вызова операции вставки. Операция удаления также объединяется с операцией поиска.

Существует большое количество модификаций предложенного подхода, ориентированных на оптимизацию поиска путем минимизации пересечений между ограничивающими кубами, улучшение процесса кластеризации данных в листовых вершинах и т. д. Наиболее известны следующие:

• R*-дерево [9]: эффективность операции поиска достигается за счет оптимизации операции вставки путем минимизации перекрытия (overlapping) между кубами узлов.

• R+^ерево [10]: процесс поиска организован эффективнее благодаря оптимизации операции вставки путем исключения overlapping. При таком подходе количество узлов в R^-дереве может быть значительно больше, а внутренняя структура сложнее. Такая структура данных занимает намного больше места в памяти, чем аналогичное R-дерево.

• X-дерево [11]: производится минимизация overlapping за счет добавления нового типа узла - суперузла, который содержит те узлы, чьи многомерные кубы пересекаются между собой. Поиск в этом узле осуществляется эффективнее, чем в обычной древовидной структуре R'''-дерева. Данный подход эффективен для многомерных пространств с размерностью больше 8.

3.2. Использование R-деревьев в задаче локализации инвариантных множеств ДС. В нашей задаче R-деревья берем в качестве представления покрытия в памяти компьютера, так как эта структура данных оперирует многомерными кубами. В отличие от задач, связанных с применением баз данных, нам достаточно хранить только координаты ячейки, поэтому потребность в хранении индексов многомерных объектов отпадает. Так как во всех реализациях алгоритма локализации инвариантных множеств рассматриваются ячейки одного размера и применяется целочисленная система координат, единица длины в которой есть размер ячейки, по нашему мнению, кажется естественным представление индексной записи в виде массива целых чисел -многомерного индекса ячейки, а листового узла - упорядоченного списка таких кортежей, элементы которого различны. Не листовой узел в нашем случае представлен как список пар: указатель на дочернюю вершину и целочисленный многомерный интервал. Если дочерняя вершина - листовой узел, тогда этот интервал равен минимальному интервалу, содержащему все номера ячеек данного листа. Иначе интервал равен минимальному интервалу, содержащему все интервалы из списка пар дочернего узла.

На практике в общем случае оценить сложность работы указанного подхода довольно затруднительно. Такая оценка напрямую связана с оценками операций поиска и добавления R-дерева. Количество операций вставки при добавлении обусловлено значениями параметров m, M и кучности ячеек в конкретном покрытии: чем больше случаев переполнений узла возникает, тем больше операций поиска и вставки инициируется. Трудоемкость операции поиска тоже зависит от покрытия: чем больше кубов узлов обладает свойством перекрытия, тем больше область для поиска.

При анализе модификаций R-дерева учитывалась не только скорость доступа к данным, но и расход памяти, требуемой для хранения ячеек покрытия в этих структурах данных. Так, в R^-дереве быстрее всех осуществляется операция поиска, так в нем нет узлов с кубами, обладающими эффектом перекрытия. Но операция вставки выполняется медленнее, чем в других подходах, так как она направлена на добавление

узлов, кубы которых не пересекаются между собой. Кроме того, Д+-дерево занимает больше компьютерной памяти, чем другие деревья, построенные на основе аналогичного покрытия. Структура Х-дерева является модификацией Д*-дерева. В предложенной задаче мы работаем с системами размерности 2 и 3. В этом случае производительность Х-дерева не выше, чем у Д*-дерева. Поэтому Д*-дерево лучше подходит при решении данной задачи.

4. Численные эксперименты. При вычислениях использовали следующую программную и аппаратную конфигурацию:

• CPU: core 2 DUO: E6420 2.13 Ghz.

• Memory 2 Gb.

• OS: Microsoft Windows XP SP 2.

Далее по тексту обозначим за ISLM + + ISLM+ с поддержкой Д*-деревьев. Количество потоков-обработчиков равно 2.

Пример 1. Рассмотрим отображение с задержкой:

xi = y,

yi = ay( 1 — x), a = 2.27.

Начальная область M = [—1.5,1.5] x [—2, 2]. Пример 2. Рассмотрим систему Икеда [12]

xi = 2 — 0.9(xcosr (x, y) — ysinr (x,y)), yi = 0.9(xcosr (x,y) — ysinr (x,y)), 6

T (x,y) = 0.4 —

1+ x2 + y

2 •

Начальная область М = [-4, 4] х [-4, 4].

Пример 3. Рассмотрим систему Дуффинга

X = y,

у = х — х3 — 0.15у.

Вычисления производились по методу Рунге-Кутта с шагом по Ь, равным 0.03. Количество итераций - 100.

В таблице приведены результаты сравнения характеристик, полученных в рассмотренных примерах.

Оравнительные характеристики работы алгортима (в мс) с оптимизациями

Система Дробление ISLM ISLM+ ISLM++

Отображение с задержкой 4 ООО ООО 16 ООО ООО 158 324 627 031 61 875 237 192 48 045 198 437

Отображение Икеда 4 ООО ООО 16 ООО ООО 203 009 815 427 135 501 540 328 107 782 429 593

Отображение Дуффинга 4 ООО ООО 16 ООО ООО Не поддержано 386 615 1417 016 350 890 1263 472

4. Заключение. В работе был представлен подход, позволивший повысить скорость обработки ДС. Такой результат достигался при помощи организации хранения

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

Литература

1. Осипенко Г. С., Ампилова Н. Б. Введение в символический анализ динамических систем. СПб.: Изд-во С.-Петерб. ун-та, 2005. 236 с.

2. Петренко Е. И. Разработка и реализация алгоритмов построения символического образа // Дифференц. уравнения и процессы управления. 2006. № 3; URL: http://www.neva.ru/journal/j/ index.html.

3. Dellnitz M., Junge O. Set oriented numerical methods for dynamical systems // Handbook of dynamical systems / ed. by B. Fiedler. 2002. Vol. 2. 270 p.

4. Ампилова Н. Б., Терентьев С. В. Применение методов интервальной арифметики к задаче построения символического образа // Дифференц. уравнения и процессы управления. 2006. № 4; URL: http://www.neva.ru/journal/j/index.html.

5. Ампилова Н. Б., Терентьев С. В. О применении интервальной арифметики при численном исследовании динамических систем // Вестн. С.-Петерб. ун-та. Сер. 10: Прикладная математика, информатика, процессы управления. 2009. Вып. 4. С. 3—11.

6. BOOST. URL: http://www.boost.org.

7. GNUPLOT. URL: http://gnuplot.sourceforge.net/.

8. Guttman A. R-Trees: A Dynamic Index Structure for Spatial Searching // Proc. ACM SIGMOD Intern. Conference on Management of Data. 1984. P. 47—57; URL: http://www.sai.msu.su/ megera/post-gres/gist/papers/gutman-rtree.pdf.

9. Beckmann N., Kriegel H. P., Schneider R. B., Seeger B. The R*-Tree: An Efficient and Robust Access Method for Points and Rectangles // SIGMOD Conference. 1990. P. 322-331; URL: http://dbs.mathematik.uni-marburg.de/publications/myPapers/1990/BKSS90.pdf.

10. Sellis T., Roussopoulos N., Faloutsos C. The R+ -Tree: a dynamic index for multi-dimensional objects. URL: http://www.vldb.org/conf/1987/P507.PDF.

11. Berchtold S., Keim D. A., Kriegei H. P. The X-tree: an Index Structure for High-Dimensional Data. URL: http://www.vldb.org/conf/1996/P028.PDF.

12. Ikeda K. Multiple-valued stationary state and its instability of the transmitted light by a ring cavity system // Opt. Commun. 1979. Vol. 30. P. 257-261.

Статья рекомендована к печати проф. Л. А. Петросяном. Статья принята к печати 24 декабря 2009 г.