Алгоритмы построения оптимальных упаковок в трехмерном евклидовом пространстве

Abstract

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

-

Постановка задачи 8i = 1; n 1; 8j = i + 1; n int B(si; r) \ int B(sj; r) = ?: Здесь int означает внутренность множества.

Сформулируем задачу о нахождении оптимальной упаковки.

Задача 1. Пусть задано ограниченное замкнутое выпуклое множество M R3 и число n 2 N. Требуется найти упаковку Un множества M , радиус r шаров в которой был бы максимальным.

Ключевым элементом решения задачи 1 является отыскание набора Sn = fsigin=1 упаковки и вычисление величины M центров шаров RM (Sn) = min min i=1;n (si; (Sn n fsig)) ; (si; @M ) ; 2 где (a; M ) расстояние от точки a до замкнутого множества M , @M граница множества M . Выражение ( 1 ) задает максимальный радиус шаров с центрами в точках из Sn, которые образуют упаковку множества M .

В предельном случае при n = 1 задача 1 сводится к построению шара максимального радиуса, вложенного в M . Заметим, что даже в этом случае ее решение требует сложных вычислений. В некотором смысле она похожа на задачу об отыскании чебышевского центра и радиуса компакта в R3 [8]. Однако в отличие от последней задачи, шар максимального радиуса, вложенный в M , может быть не единственным, даже если M выпуклый многогранник. Обозначим (M ) = fm 2 M : (m ; @M ) = max (m; @M )g множеm2M ство всех точек множества, в которых достигается максимум расстояния до его границы. По построению, любая точка из (M ) есть центр шара наилучшей 1-упаковки в M .

Для удобства при построении алгоритмов решения задач 1 формулу ( 1 ) можно записать в виде где

RM (Sn) = min 'i(si);

i=1;n 'i(x) = min 1 2 j=1;n(j6=i) kx sik; (x; @M ) ; i = 1; n; min ( 1 ) ( 2 ) ( 3 ) здесь k k означает норму вектора.

Соответственно, процесс отыскания множества Sn центров шаров упаковки можно свести к поэтапной максимизации функций 'i(x), i = 1; n; на множестве M . При этом i меняется от 1 до n, точки Sn n fsig считаются фиксированными, а si строится как точка локального максимума функции ( 3 ) при заданном i.

К форме задачи 1 можно привести различные постановки, известные в геометрии. В частности, И. Ньютон сформулировал задачу о так называемом контактном числе шаров максимальном количестве шаров единичного радиуса, которые могут одновременно касаться одного такого же шара в R3 (при условии, что шары пересекаются только в граничных точках). Можно разместить 12 шаров единичного радиуса с центрами в вершинах правильного икосаэдра так, чтобы они касались одного шара. Но у Ньютона не было строгого доказательства, что нельзя разместить 13 шаров так, чтобы они имели точку пересечения с некоторым одним шаром. В терминах данной статьи эту задачу можно сформулировать так. Задано множество M = fm : 1 6 kmk 6 3g. Требуется решить задачу о построении оптимальной упаковки при n = 13 и вычислить радиус ее шаров r . В работе [9] доказана оценка r < 1, которая означает, что контактное число равно 12, а не 13. При этом при n = 14 радиус шаров равен r 0:98, то есть лишь ненамного меньше 1 [10].

Методы построения наилучших упаковок Применительно к ограниченным множествам оптимальные упаковки можно строить аналитически только для самых простых случаев. Одним из них является задача 1 при n = 1 для выпуклого многогранника M . Как уже было отмечено выше, в этом случае требуется отыскать произвольную точку из множества (M ). Из свойств выпуклых многогранников известно, что шар B(x; r) наибольшего радиуса, вложенный в M , расположен так, что x 2 co(@B(x; r) \ @M ), где co M означает выпуклую оболочку множества M . Действительно, если бы данное условие не выполнялось, то можно было бы сдвинуть точку x так, чтобы расстояние от нее до ближайшей к ней грани многогранника увеличилось, и получить точку x , являющуюся центром шара большего радиуса, вписанного в M . Соответственно, можно искать точку x 2 M как точку, равноудаленную от некоторых четырех граней многогранника. В свою очередь, эту задачу можно свести к решению системы линейных уравнений, задающих плоскости, равноудаленные от плоскостей, содержащих грани M . Если найденная при их решении точка x лежит в M , ее проекции на плоскости принадлежат соответствующим граням и x не имеет более близких проекций на остальные грани, то она принадлежит (M ). Данный метод считать некоторым аналогом метода перебора симплексов из вершин многогранника при отыскании его чебышевского центра [11]. Условие выпуклости является существенным, поскольку для произвольного невыпуклого многогранника шар наибольшего радиуса, вложенный в него, может пересекать своей границей ребра или вершины многогранника.

Заметим, что плоскость , равноудаленная от двух плоскостей i и j, содержащих некоторые две грани i и j, является единственной только в том случае, если i и j параллельны. В противном случае найдутся две ортогональные бисекторные плоскости и , проходящие через прямую i \ j, точки которых лежат на равном расстоянии от i и j. Однако поскольку по условию M выпуклый многогранник, то он полностью лежит в одном из полупространств, ограниченных плоскостью, проходящей через его грань. Следовательно, всегда можно указать ту единственную плоскость из пары ; , непустое подмножество точек которой вложено в int M . В случае рассмотрения смежных граней обе бисекторные плоскости имеют непустое пересечение с M , но для одной из них это пересечение совпадает с одной из граней M , поэтому его точки не могут являться центрами шаров B(r; s) M радиуса r > 0. Решение задачи 1 реализовано на базе методов вычислительной геометрии Определение 2. Пусть задано компактное множество M и множество Sn, состоящее из n точек. Областью Дирихле [12, c. 305] точки si 2 Sn во множестве M называется подмножество

Di(M; Sn) = fm 2 M : km sik = (m; Sn)g: В случае если Sn M , все области Дирихле Di(M; Sn); i = 1; n, непустые компакты. Если M выпуклое множество, то и области Дирихле в нем выпуклые при произвольном Sn. Если M выпуклый многогранник, то множества Di(M; Sn); i = 1; n, тоже являются по построению выпуклыми многогранниками. Для построения областей Дирихле авторы используют так называемую трехмерную диаграмму Вороного множества Sn, то есть геометрическое место точек, у которых имеются два или более ближайших элемента из Sn (подробнее см. [13]). Диаграмма Вороного в общем случае разбивает пространство на области, лежащие ближе к одной из точек si, нежели к sj; j = 1; n; j 6= i:

Лемма 1. Пусть задано выпуклое множество M и множества Sn M из n попарно различных точек. Тогда справедливо 8i = 1; n 'i(si) = ( 4 ) ( 5 ) Доказательство. Без ограничения общности полагаем i = 1. Докажем неравенство Допустим, оно не выполняется. Тогда найдется точка m 2 @Di(M; Sn) такая, что 'i(s1) 6

s1; @D1(M; Sn) : ks1

mk < 'i(s1): По построению области Дирихле, точка m может либо принадлежать @M , либо лежать на равном расстоянии от s1 и некоторой другой точки из Sn. Первый случай невозможен, поскольку тогда значение величины (s1; @M ) не может быть больше, чем ks1 mk, а значит, и значении функции 'i(s1) тоже. Если реализовался второй случай, то можно указать точку из Sn (без ограничения общности полагаем, что s2) такую, что ks1 mk = ks2 mk: Тогда в силу неравенства треугольника имеет место ks1 s2k 6 ks1 mk + ks2 mk = 2ks1 mk; следовательно, 'i(s1) 6 ks1 s2k=2 6 ks1 Докажем неравенство

mk. Получилось противоречие. Допустим, оно не выполняется. Тогда найдется точка m , либо принадлежащая @M , либо являющаяся серединой отрезка, с концами соответственно в s1 и в некоторой другой точке из Sn такая, что ( 6 ) 'i(s1) > Первый вариант невозможен, поскольку, если точка m лежит на границе множества M и не принадлежит @D1(M; Sn), то она лежит в некоторой другой области Дирихле Dj(M; Sn); j > 1. Следовательно, можно построить отрезок [s1; m ], у которого конец m лежит вне D1(M; Sn). Тогда найдется точка m 2 [s1; m ] \ @D1(M; Sn). В итоге имеет место неравенство Если m не принадлежит @D1(M; Sn) , то она, как и в предыдущем случае, лежит в некоторой области Дирихле Dj(M; Sn); j > 1, такой, что ks1 m k > ksj m k. Соответственно, отрезок [s1; m ] минимум в одной точке m пересекает множество @D1(M; Sn). Получается противоречие, как и в предыдущем случае. Из оценок ( 5 ) и ( 6 ) следуют равенства ( 4 ).

На базе леммы 1 разработан алгоритм итерационного улучшения множества точек Sn M с позиции максимизации значения величины ( 2 ), представляющий из себя переработку алгоритма из работы [14] для случая трехмерного пространства.

Алгоритм 1. 1. Переменная i циклически меняется от 1 до n. 2. Строится область Дирихле Di(M; Sn) с помощью разбиения многогранника M на части диаграммой Вороного. 3. Посредством перебора всех четверок из числа граней многогранника Di(M; Sn) находится точка bsi 2 (Di(M; Sn)) центр шара максимального радиуса, вложенного в Di(M; Sn). В случае если (Di(M; Sn)) содержит более, чем 1 элемент, то bsi выбирается из них с помощью генератора случайных чисел. 'j(sj) = ksi 2

sjk : 'i(si) > 'i(si): 8j 2 J 8j 2 J 'j(sj) > 'i(si): 'j(sj) > 'i(si): 8j 2 1; n 'j(sj) > 'j(sj); Заметим, что неравенство 'j(sj) > 'j(sj) может иметь место только при j 2 J .

Действительно, значение функции 'i(si) не может быть меньше, чем 'i(si) в силу леммы 1. Значит Если же j 6= i и j 2= J , то для точки sj найдется либо точка на границе, лежащая к ней ближе, чем ksi sjk=2, либо точка sk 2 Sn(k 6= i; k 6= j) такая, что ksk sjk < ksi sjk. Значит, можно записать оценку Из определения функций 'j следует оценка 8j = 1; n (j 6= i; j 2= J ) 'j(si) > 'j(si): 8j = 1; n (j 6= i) 'i(si) 6 ksi sjk=2: По построению множества J имеем С учетом неравенства ( 9 ) запишем Из ( 9 ), ( 10 ) и ( 11 ) вытекает оценка Обозначим через J множество номеров тех точек из Sn, для которых выполняется равенство ( 7 ) ( 8 ) ( 9 ) ( 10 ) ( 11 ) равномерное распределение на отрезка [ "; "], причем " существенно меньше l. Для найденной точки si выполняется проверка условия, принадлежит ли она M , и не совпадает ли с каким-либо уже записанным элементом Sn. При стохастической генерации массива точек возможна такая их конфигурация, которая далека от оптимальной, однако локально устойчива относительно применения алгоритма 1. Например, если M куб, а точки si 2 Sn0 расположены на отрезке, соединяющем центры противоположных граней. Однако на практике программа запускается многократно при различных значениях Sn0, что позволяет отбрасывать подобные упаковки с малым радиусом шаров RM (Sn). 4

Примеры Авторами разработан в пакете MATLAB программный комплекс, реализующий построение аппроксимаций оптимальных упаковок для выпуклых многогранников. В нем использованы процедуры вычислительной геометрии, которые раньше применялись для решения прикладных задач, в частности, построения функции евклидового расстояния до множества, сингулярных линий в задачах быстродействия [15] и оптимальных покрытий [11].

Пример 1. Требуется построить наилучшие упаковки Un для куба M = f(x; y; z) 2 R3 : jxj 6 1; jyj 6 1;jzj 6 1g при числе шаров n, равном 20 и 30.

У куба 6 граней i; i = 1; 6, являющиеся квадратами. При генерации первоначального приближения множества Sn применялись стохастические процедуры с ограничением на координаты точек si; i = 1; n, чтобы они не превышали по модулю 0:9.

При n = 20 радиус шаров аппроксимации оптимальной упаковки в куб M равен r 0:3412. Ребра куба M и упаковка шаров U20 представлены на рис. 1 и 2. Массив центров шаров

S20

f( 0:5312; 0:6547; 0:6300); ( 0:6554; 0:5032; 0:0334); ( 0:4099; 0:6584; 0:6584); ( 0:6551; 0:0239; 0:6551); ( 0:2232; 0:0387; 0:0078); ( 0:6587; 0:0223; 0:6587); ( 0:3128; 0:6558; 0:6558); ( 0:6566; 0:5764; 0:0656); ( 0:6588; 0:6588; 0:6141); (0:1595; 0:6557; 0:6557); (0:0163; 0:6556; 0:0179); (0:2751; 0:6566; 0:6566); (0:0336; 0:0602; 0:6553); (0:4485; 0:2764; 0:0266); (0:0378; 0:0009; 0:6514); (0:6312; 0:6312; 0:6312); (0:0114; 0:6465; 0:3916); (0:6392; 0:4346; 0:6392);

(0:6554; 0:2000; 0:5110); (0:6546; 0:3721; 0:1572)g: При n = 30 радиус шаров аппроксимации оптимальной упаковки равен r упаковка шаров U30 представлены на рис. 3 и 4. Массив центров шаров 0:2953. Ребра куба M и S30

f( 0:6895; 0:6895; 0:5598); ( 0:6992; 0:6992; 0:2497); ( 0:3158; 0:5792; 0:6984); ( 0:1844; 0:3548; 0:6991); ( 0:5385; 0:2035; 0:2376); ( 0:2228; 0:1535; 0:2774); ( 0:6981; 0:1127; 0:6782); ( 0:6964; 0:1627; 0:6964); ( 0:2388; 0:3202; 0:2540); ( 0:7000; 0:2081; 0:1682); ( 0:1300; 0:2759; 0:6966); ( 0:4145; 0:7050; 0:6991); ( 0:7009; 0:7009; 0:1726); ( 0:3071; 0:6978; 0:2850); ( 0:6977; 0:4912; 0:6977); (0:3124; 0:6974; 0:6974); (0:5108; 0:6970; 0:1253); (0:1874; 0:6977; 0:3864); Пример 2. Требуется построить наилучшие упаковки Un для октаэдра M = f(x; y; z) 2 R3 : jxj + jyj+ +jzj 6 1g при числе шаров n, равном 20 и 30.

У октаэдра 8 граней i; i = 1; 8, являющихся правильными треугольниками. Октаэдр дуален кубу, то есть вершины одного из этих многогранников центры граней другого. При генерации первоначального приближения множества Sn применялись стохастические процедуры с ограничением на координаты точек si; i = 1; n, чтобы суммы модулей абсциссы, ординаты и аппликаты не превышали 0:9. Рис. 1: Аппроксимация U20 наилучшей упаковки ку- Рис. 2: Аппроксимация U20 наилучшей упаковки куба M 20-ю шарами: вид 1. ба M 20-ю шарами: вид 2. Рис. 3: Аппроксимация U30 наилучшей упаковки куба M 30-ю шарами: вид 1. Рис. 4: Аппроксимация U30 наилучшей упаковки куба M 30-ю шарами: вид 2. При n = 30 радиус шаров аппроксимации оптимальной упаковки равен r упаковка шаров U30 представлены на рис. 7 и 8. Массив центров шаров 0:1570. Ребра куба M и ( 0:3074; 0:2362; 0:1956); ( 0:3363; 0:3554; 0:0326); ( 0:2296; 0:4976; 0:0129);

( 0:0971; 0:1007; 0:0055); (0:5502; 0:1782; 0:0003); (0:0718; 0:0609; 0:2568); ( 0:3080; 0:1648; 0:2488); ( 0:1042; 0:4095; 0:2089); ( 0:2363; 0:1368; 0:3436); ( 0:0369; 0:2307; 0:0036); (0:5927; 0:1337; 0:0000); (0:2537; 0:0729; 0:0041); (0:1490; 0:1599; 0:2269); (0:0000; 0:0000; 0:7171); (0:0000; 0:3663; 0:3502); ( 0:0853; 0:6525; 0:0003); (0:4408; 0:0008; 0:2681); (0:1666; 0:1004; 0:4406);

(0:2026; 0:5228; 0:0000); (0:0894; 0:5580; 0:0924); (0:2259; 0:3656; 0:1337)g: Рис. 5: Аппроксимация U20 наилучшей упаковки октаэдра M 20-ю шарами: вид 1. Рис. 6: Аппроксимация U20 наилучшей упаковки октаэдра M 20-ю шарами: вид 2. Рис. 7: Аппроксимация U30 наилучшей упаковки октаэдра M 30-ю шарами: вид 1. Рис. 8: Аппроксимация U30 наилучшей упаковки октаэдра M 30-ю шарами: вид 2.

В процессе построения упаковок в каждом примере параметры окончания работы программного комплекса были равны r = 10 4; h = 10 3. Для каждого из значений n было выполнено от 10 до 15 запусков с различными начальными значениями Sn0. Общее число циклов алгоритма 1 варьировалось от 39 до 77. Во всех случаях n-сеть пришла к установившемуся значению. При графическом представлении результатов упаковки показаны под разным углом: в одном случае фронтальная плоскость близка к плоскости xOz, в другом к yOz. Авторами разработан программный комплекс построения оптимальных упаковок шаров в многогранники. Он основан на геометрических алгоритмах максимизации функции ( 1 ), которые итерационно изменяют заданное множество точек Sn. Для построения начального приближения массива центров шаров Sn использованы стохастические процедуры, учитывающие геометрию многогранника. Показано, что алгоритмы являются монотонными, то есть не уменьшают радиус шаров. Проведено численное моделирование примеров; построены аппроксимации наилучших упаковок шаров в куб и в правильный октаэдр. Выполнена визуализация результатов. 6

Благодарности Список литературы

Работа была выполнена при финансовой поддержке РФФИ (проект №16-31-00356-мол_а) и Программы фундаментальных исследований УрО РАН, проект №15-16-1-13.

Algorithms of optimal packing construction in a 3-dimensional Euclidian space

Pavel D. Lebedev, Alexandr A. Uspenskii Krasovskii Institute of Mathematics and Mechanics (Yekaterinburg, Russia) Keywords: optimal packing, polyhedron, Dirichlet domain, distance function, stochastic algorithm.

The problem of optimal equal balls packing construction into the compact set M in a 3-dimensional Euclidian space is studied. Criteria of optimality is considered radii of the balls. Numerical methods based on splitting M into subsets and calculating balls inscribed into them are designed. The theorem about suggested algorithm properties is proved. Constructing optimal packing approximations is produced with large quantity of its elements. The results are presented as 3-dimensional illustrations.

References

[1] Krasovskii N. N. , Subbotin A. I. Positional Differential Games. Nauka , Moscow. 1974 (in Russian). = Красовский Н. Н. , Субботин А . И. Позиционные дифференциальные игры . М.: Наука , 1974 .

[2] Kazakov A. L. , Lempert A. A. An Approach to Optimization in Transport Logistics . Autom. Remote Control . 72 ( 7 ): 1396 - 1404 , 2011 .

[3] Yamada S. , Kanno J. , and Miyauchi M. Multi-sized sphere packing in containers: optimization formula for obtaining the highest density with two different sized spheres . IPSJ Online Transactions . 4 : 126 - 133 , 2011 .

[4] Stoyan Y. , Yaskov G . Packing congruent spheres into a multi-connected polyhedral domain Intl . Trans. in Op. Res . 20 : 79 - 99 , 2013 . DOI: 10 .1111/j.1475- 3995 . 2012 . 00859 .x

[5] Conway J. H. , Sloane N. J. A. Sphere Packing , Lattices and Groups Springer, New York, 1988 .

[6] To¨th L. F. Lagerungen in der Ebene, auf der Kugel und im Raum Springer. Springer, Berlin, 1957 .

[7] Sloane N. J. A. Scientific American . 250 ( 1 ): 116 - 125 , January 1984 . DOI: 10 .1038/scientificamerican0184- 116

[8] Garkavi A. L. On the Chebyshev center and convex hull of a set . Uspekhi Mat. Nauk . 19 ( 6 ), 139 - 145 , 1964 .

[9] Schu¨tte K., van der Waerden B. L. Das Problem der dreizehn Kugeln Math . Ann. 125 ( 1 ): 325 - 334 , 1953 . DOI: 10 .1007/BF01343127.

[10] Musin O. R. , Tarasov A. S. The Tammes problem for N = 14 . E-print, 2014 , arXiv: 1410 .2536 [math.MG]

[11] Lebedev P. D. and Ushakov V. N. Algorithms for the construction of an optimal cover for sets in threedimensional Euclidean space . Tr. IMM UrO RAN . 21 ( 2 ): 276 - 288 . 2015 .

[12] Brusov V. S. , Piyavskii S. L. A computational algorithm for optimally covering a plane region . USSR Comp. Math. Math. Phys . 11 ( 2 ): 17 - 27 , 1971 .

[13] Mestetsky L. M. Continuous Morphology of Binary Images: Shapes, Frames, Circulars . Fizmatlit, Moscow, 2009 (in Russian). = Местецкий Л . М. Непрерывная морфология бинарных изображений: фигуры, ске- леты, циркуляры . М.: Физматлит , 2009 .

[14] Kazakov A. L. and Lebedev P. D. Algorithms of Optimal Packing Construction for Planar Compact Sets . Numerical Methods and Programming . (Vychislitel'nye Metody i Programmirovanie) . 16 ( 3 ): 307 - 317 , 2015 (in Russian). = Казаков А. Л. , Лебедев П . Д. Алгоритмы построения оптимальных упаковок для ком- пактных множеств на плоскости . Вычислительные методы и программирование . 16 ( 3 ): 307 - 317 , 2015 .

[15] Lebedev P. D. , Uspenskii A. A. On the set of limit values of local diffeomorphisms in wavefront evolution . Proceedings of the Steklov Institute of Mathematics . 272 ( 1 ): 255 - 270 , 2011 .