Математическая модель и метод решения задачи упаковки гомотетических одинаково ориентированных эллипсоидов

УДК 519.85

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ И МЕТОД РЕШЕНИЯ ЗАДАЧИ УПАКОВКИ ГОМОТЕТИЧЕСКИХ ОДИНАКОВО ОРИЕНТИРОВАННЫХ ЭЛЛИПСОИДОВ

ХЛУД О.М.,СУББОТА И.А.РОМАНОВА Т.Е.

Рассматривается задача упаковки гомотетичных одинаково ориентированных эллипсоидов в контейнере минимального объема. В качестве контейнера выбирается прямоугольный параллелепипед или эллипсоид. Строится математическая модель в виде задачи нелинейного программирования с использованием phi-функций. Предлагается эффективный алгоритм решения, использующий гомотетические преобразования эллипсоидов и оптимизационную процедуру, котрая позволяет сократить вычислительные ресурсы. Приводятся результаты численных экспериментов.

Ключевые слова - упаковка, гомотетичные эллипсоиды, phi-функции, математическая модель, нелинейная оптимизация

Введение

Задачи упаковки и раскроя (Cutting and Packing) относятся к классу NP-сложных и являются предметом исследования вычислительной геометрии, а методы их решения - новым направлением теории исследования операций.

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

Многие публикации посвящены решению данного класса задач (см., например, [1-4]).

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

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

Имеется набор гомотетичных одинаково ориентированных эллипсоидов Е;, 1 е IN = {1,2,..., №}, заданных полуосями а;, Ь;, С; в собственной системе координат O¡xyz . Полагаем, что начало собственной системы координат Е; находится в центре его симметрии. Положение Е; в пространстве R3 определяется переменным вектором трансляции 3 3

V; = (х^у^^ е где R - евклидово пространство.

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

т.е. р = (А, В, С), либо эллипсоид е , заданных полуосями А, В, С в фиксированной системе координат Oxyz с переменным коэффициентом гомотетии X, т.е. р = (X).

В зависимости от вида контейнера рассматриваются следующие функции цели: F = А • В • С, если О = Р; Б = X, если О = Е.

Задача. Упаковать набор одинаково ориентированных гомотетичных эллипсоидов {Е;(и;), 1 е в контейнер О так, чтобы функция цели б достигала своего минимального значения.

2. Математическая модель

Основными ограничениями поставленной задачи являются:

- непересечение эллипсоидов, т.е.

;)пintЕj

- включение эллипсоида в область размещения, т.е.

Е^ 1) сОо ЫЕ^ 1) п Ш О * =0, 1 е 1М, (2) где О * = R3 \Щ О.

Как известно [5], конструктивным средством математического моделирования ограничений (1), (2) является метод рЫ-функций Стояна [6]. В терминах рЫ-

функций ограничение (1) имеет вид ф^(Л1 ьvj)^0,

где

Ч (Xj - Xi)2 (yj - yi)2 (Zj - Zi)2 ,

Ф ij

j j

- phi-функция для эллипсоидов Ej(v ¡) и Ej(v j), а ограничение (2) описывается неравенством

Ф¡(V ¡) > 0 , где Ф¡(V ¡) - рЫ-функция для эллипсоида Ej(v¡)и объекта . При этом:

если О = Р , то

Фi(vi) = min хk(vi) k=1,...,6 :

где

Xi(v i) = - xi + A - аь х 2(v i) = -yi + B - bh

X3(vi) = -zi + C-cb X4 (vi) = xi + A -ab

х 5(v i) = yi + B - bi, X6(v i) = zi + C - ci;

если q = e , то

Oi

(XA - a;)2 (XB - b;)2 (XC - 0;)2 ' W

при условии, что Ei и Е гомотетичны.

Математическая модель поставленной задачи может быть представлена так:

min F

W = {u e R ст : Ф ij

где u = (p, v i,..., v n ) - вектор переменных,

p = (A, B, C), если Q = P , p = (X), если Q = E, R ст -арифметическое евклидовое пространство размерности ст, Z > 0 - система дополнительных ограничений на метрические характеристики контейнера Q .

Задача (1)-(2) является многоэкстремальной задачей нелинейного программирования, где целевая функция линейна или квадратична, а множество допустимых решений задается с использованием квадратичных функций (3), (5) и кусочно-линейных функций

(4).

3. Метод решения

Для поиска локальных минимумов задачи (6)-(7) используется подход, в основе которого - метод мультистарта и оптимизационная процедур а, включающая поиск допустимых стартовых точек и локальную оптимизацию, которая является развитием алгоритмов, предложенных в статьях [7, 8]. В качестве локально-оптимального решения выбирается наилучший из полученных локальных экстремумов.

3.1. Пошаговый алгоритм

Полагаем, что стартовые размеры контейнера Q0 достаточно большие и гарантируют размещение в нем эллипсоидов Ei. Для контейнера Q = P в качестве стартовых параметров выбираем A0, B0, C0 , а для Q = E выбираем X 0. При этом, не теряя общности, полагаем, что ai >a2 >.... > aN-i > aN.

Пусть коэффициенты гомотетии X; эллипсоидов Е;, i е I м, являются переменными. Тогда и = (V, X) е R ° - вектор переменных, где X = (X!,..., Xм) е Rм .

Шаг 1. Выбираем точку и0(1) = (V0,X0), где V = (х?,у?^?) еП0 выбирается случайно, X 0 = (X!,..., X N), X 0 = 0, i е I м . В качестве стартовой точки для дальнейшей оптимизации выбираем и0(1) и переходим к шагу 2.

Шаг 2. Решаем задачу

N

max £ ai X i ueW'cR ст i=1

W' = {ue Rст :Фij(u) >0,Ф;(u) >0,

0¿Xi < 1, j>¡, iе 1К}.

В результате ее решения получаем точку глобального максимума и0(2) = ^0(2), X0(2)). В качестве стартовой точки для дальнейшей оптимизации выбираем и0(2) и переходим к шагу 3.

Шаг 3. Решаем задачу (5)-(6). В результате решения получаем точку локального минимума

и0(3) = (p0(3),v0(3)).

Шаг 4. В качестве стартовой точки для дальнейшей оптимизации выбираем точку и0(3) = ("у0(3), X 0(3)) при фиксированных размерах контейнера. Решаем вспо-

могательную задачу вида

N

max £

ueW'tR ст i=1

W" = {u e Rст : Фij(u) > 0, Ф;(u) > 0, a- <aiXi <a+, i< j = 1,...,N}, a- = min{ai,i e IN} , a + = max{ai, i e IN}.

В результате решения задачи (10)-(11) получаем точку локального максимума и 0(4) = (-у0(4), X 0(4)) .

Шаг 5. Ранжируем по убыванию аi = X, i е IN . Формируем последовательность (^^^...Л^Лм), такую, что а^ >а\2 >....>а¡м-1 >а.

Осуществляем сравнение a¡ и а^ : если X0(4) > 1, то

X 0(5) = 1 1 0(4) ,

полагаем -1, если X- < 1, то полагаем

^ ' i '

X 0(5) =Х0(4). Формируем точку и0(5) = (у°(5), X 0(5)) и переходим к следующему шагу.

Шаг 6. Стартуя из точки и0(5) = (у0(5), X0(5)), решаем задачу (8)-(9). В результате решения получаем точку глобального максимума и0(6) = (у0(6), 6 0(6), X 0(6)). В противном случае переходим к шагу 1.

Стартуя из точки и0(6), решаем задачу (6)-(7) и получаем точку локального минимума

и'(6) = (р'(6),у'(6)).

Шаг 7. Генерируем случайно точку q1. Если

ql йЕ1, 1 е 1М, q1 еО, то в качестве параметров размещения дополнительного эллипсоида с полуосями аq1 = а 1, ьql = Ь1 , сql = с1 и коэффициентом

гомотетии Xql = е выбираем точку ql.

Стартуя из точки и0(7) = (у'(6), V0 ^, X 0 ^), при фиксированных размерах контейнера решаем следующую вспомогательную задачу:

тах а q1 X q1

uеWсRа 41

Wx= {иеR° :х

где и = (V, V ql, X ql) е R

3.2. Алгоритм локальной оптимизации c преобразованием области допустимых решений

В основе алгоритма лежит оптимизационная процедура LOFRT, предложенная в [9] для задачи оптимального размещения эллипсов, и состоит в следующем.

Пусть и(0) е W - допустимая точка. Около каждого эллипсоида Е; описывается сфера 81 радиуса а;, 1 е IN . Для каждой сферы 8; строится «индивидуальный» кубический контейнер 01 з Si з Е1, длина стороны которого составляет 2 • (а 1 + е), 1 е IN . При этом 81, Е1 и 01 имеют один и тот же центр (х10,у0^0), а стороны 01 параллельны соответствующим сторо-

N

нам О . Выбирается е = ^ Ь1 / N. Далее фиксируется 1=1

положение каждого индивидуального контейнера О1 (рис.1). Каждый эллипсоид Е1 может размещаться только внутри индивидуального контейнера О1 .

/ £

V = (V1,..., V N ),

Х(и) = тт{фх(и), Ф1 (и), т = 1,...,т,1 е т = 0.5• N• ^ +1) .

В результате решения задачи (12)-(13) получаем точ-

0(7) _ , 0(7) 0(7) л 0(7) ч ку локального максимума и = (у , ,Xq1 ).

Осуществляем поиск ^ aq1 , ^ ЬЯ1 ,

] е IN . Если такой эллипсоид Ej существует, то

принимаем ^ = v q1 по аналогии с алгоритмом, предложенным в работе [7].

Формируем точку: и'(7) = (р'(6),у'(7)). Стартуя из

точки и'(7) = (р'(6),у'(7)), решаем задачу (6)-(7). Получаем точку локального минимума и''(7) = (р'' (7),у''(7)), причем и* = и"(7) = (р*,у*, 6*).

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

Рис. 1. Формирование индивидуального контейнера

Далее определяется система дополнительных огр ани-чений на вектор трансляции У1 для каждого эллипсо-

ида Е1 в виде Ф81 °и > 0 , 1 е I , где

Ф81°н = тш{—х1 +х? +е, —у1 + У00 +е, -zi +z11 +е, Х1 — х0 +е, у1 —у0 +е, Zi —z° +е}

- рЫ-функция для 81 и О* = R3 \int О11.

Заметим, что неравенство ф 8 > 0 эквивалентно системе линейных неравенств — х1 +х0 +е> 0,

—у1 +у0 +е >0 , —zi +z° +е , х1 — х0 +е>0, у 1 —у0 +е> 0, zi —z° +е.

Очевидно, что если контейнеры 01 и О j не имеют общих внутренних точек, т.е. Ф°1ОJ > 0, тогда не следует накладывать ограничение на непересечение эллипсоидов Е1 и Ej. Например, для эллипсоидов Е1 и Е7, Е4 и Е8, Е1 и Е8 (рис. 2).

I I

С 1 —-V4 >

у гТ^

ttt^vIO

С?"

L4-

а б

Рис. 2. Размещение эллипсоидов Е; и их индивидуальных контейнеров О; : а - для О = Р , б - О = Е

Далее формируется область допустимых решений

Wi ={u eR0 : Ф^ > 0,(i,j) eSb Фi > 0,i eS2, Ф^"1* >0,ielN,l>l0 -e,w > w0 -e,h>h0 -e}:

min

uwk eWkcR<

F(Uwk),

41 = {(i,j): Ф^"4 < 0, i > j e IN},

ik2 = {i:Ф°ki° <0,ieIn}-

= {(i, j):ф"1" <0}, E2 = {i:ф""' <0}.

Другими словами, из системы, которая описывает область W, исключаем phi-неравенства для тех пар эллипсоидов, у которых индивидуальные контейнеры не пересекаются. При этом добавляем вспомогательные неравенства фsi"h > о , описывающие условие включения S; в соответствующий индивидуальный контейнер "j, i е IN .

Затем осуществляем поиск точки локального минимума uWl для подзадачи вида

min F

uw1eWi cR a .

Точка uWl является стартовой точкой u(1) на второй

итерации оптимизационной процедуры. На данном этапе вновь определяются все пары эллипсоидов с непересекающимися индивидуальными контейнерами, формируется соответствующая подобласть W2

(аналогично Wi) и вычисляется точка локального

*

минимума uW2 е W2, которая используется в каче-

стве стартовой точки uv у на третьей итерации, и т.д.

В общем случае, на k-м шаге, стартуя из точки u(k 1), решается подзадача вида

Wk ={u eR 0 : фк > 0, (i,j) eSki, Фк > 0,i eSk2, Ф^"* > 0

T11* * 11*

i eIN, 1 >lwk -6, w >wwk -6, h > hwk -sb

SiQki

Итерационная процедура заканчивается, когда

F(uWk)=F(uWk+1).

* (k)* * о

Точка u = uw = uwk e R является точкой локаль-

* о

ного минимума задачи (1)-(2), где uWk e R является точкой локального минимума на последней итерации.

Таким образом, для 0(n2) пар эллипсоидов в контейнере LOFRT процедура позволяет осуществлять проверку phi-функций, в общем случае, только для

O(n) пар эллипсоидов (это зависит от размеров эллипсоидов и величины е).

Параметр е обеспечивает баланс между количеством неравенств в каждой подзадаче нелинейного программирования (14)-(15) и числом подзадач, которые необходимо решить для поиска локально-оптимального решения задачи (6)-(7).

Таким образом, предложенный алгоритм позволяет свести задачу (6)-(7) с количеством неравенств 0(n2) к последовательности задач с количеством неравенств O(n). Это приводит к значительному сокращению вычислительных ресурсов при решении задач нелинейного программирования.

4. Численные эксперименты

Чтобы продемонстрировать эффективность предлагаемого подхода, приведем ряд примеров с использованием описанного выше алгоритма. Эксперименты проводились на компьютере Intel(R) Core(TM) i7-3630QM. Для поиска локальных минимумов использовался солвер FindArgMin пакета Wolfram Mathematica 9.

Во всех примерах осуществляется поиск 10 локальных минимумов.

В первых двух примерах рассматривается N=20 эллипсоидов, имеющих следующие размеры:

(aj, b!, c 1) = (30, 10, 10) ,(a2,b2,c2) = (9, 3, 3) , (Яз,Ъз,Сз) = (7.5, 2.5, 2.5) ,{(ai,bi,ci) =(6, 2, 2),

i = 4,...,11},{(ai,bi,ci) = (3, 1, 1) ,i = 12,..., 20} .

Пример 1. Размещение эллипсоидов в контейнере-

параллелепипеде, соответствующее точке локального

*

минимума u , приведено на рис. 3. Контейнер имеет объем F(u*) = 3213.92.

Использование LOFRT процедуры позволило значительно сократить время решения задачи. В частности, для примера 3 время решения без использования алгоритма - 13 часов, с его применением - 5 часов.

Пример 4. Размещение эллипсоидов в контейнере-

эллипсоиде, соответствующее точке локального ми*

нимума и , приведено на рис. 6. Контейнер имеет коэффициент гомотетии F(u*) = 0,45775.

Рис. 3. Локально-оптимальное размещение эллипсоидов для примера 1

Пример 2. Размещение эллипсоидов в контейнере-

эллипсоиде, соответствующее точке локального ми*

нимума и , приведено на рис. 4. Контейнер имеет коэффициент гомотетии F(u*) = 0,32099.

Рис. 4. Локально-оптимальное размещение эллипсоидов для примера 2

Далее рассматривается упаковка N=50 эллипсоидов, имеющих следующие размеры:

(а1,Ь1,с1) = (60, 20, 20) ,(а2,Ь2,с2) = (30,10,10),

(а3,Ь3,с3) = (18,6,6) ,{(а1,Ь1,с1) = (10.5, 3.5, 3.5) , 1 = 4,...,8},{(а1,Ь1,С1) = ( 9, 3, 3), 1 = 8,...,13},

{(а1,Ь1,с1) = (8.25, 2.75, 2.75), 1 = 14,...,17} , {(а1,Ь1,с1) = (7.5,2.5,2.5),1 = 18,...,24} ,

{(а1,Ь1,с1) = =(6,2,2),1 = 25,...,40} ,

{(а1, Ь1,С1) = (3,1,1) ,1 = 41,..., 50} .

Пример 3. Размещение эллипсоидов в контейнере

О = Р , соответствующее точке локального миниму-

*

ма и , приведено на рис. 5. Контейнер имеет объем

*

F(u ) = 33 874.5 .

Рис. 5. Локально-оптимальное размещение эллипсоидов для примера 3

Рис. 6. Локально-оптимальное размещение эллипсоидов для примера 4

Пример 5. Рассматривается упаковка N=50 эллипсоидов, имеющих следующие размеры:

{(а1,Ь1,с1) = (30,10,10) ,1 = 1,2} , (а3,Ь3,с3) = (18,6,6), {(а1,Ь1,с1) = (10.5,3.5,3.5),1 = 4,...,8} , {(а1,Ь1,с1) =

= ( 9, 3, 3), 1 = 8,...,13},{(аь Ьь с1) = = (8.25, 2.75, 2.75) , 1 = 14,...,17},{(а1,Ь1,с1) =

= (7.5,2.5,2.5), 1 = 18,..., 24} , {(аьЬьс1) = (6,2,2),

1 = 25,..., 40} ,{(а1,Ь1,с1) = (3,1,1), 1 = 41,...,50}.

Размещение эллипсоидов в О = Р , соответствующее точке локального минимума и* , приведено на рис. 7. Контейнер имеет объем F(u*) = 8030.25 .

Рис. 7. Локально-оптимальное размещение эллипсоидов для примера 5

Использование алгоритма LOFRT позволило уменьшить время решения задачи для примера 5 на 10 часов (время решения без использования алгоритма - 18 часов, с его применением - 8 часов).

Пример 6. Рассматривается упаковка 75 гомотетичных одинаково ориентированных эллипсоидов, имеющих следующие размеры:

{(а1,Ь1,с1) = (10.5, 3.5, 3.5), 1 = 1, ...,15},

{(а1,Ь1,с1) = ( 9,3,3) ,1 = 16,...,30},

{(ai,bi,ci) = =( 7.5,2.5,2.5) ,i = 31,...,45}, {(ai,bi,ci) =(6,2,2), i = 46, ...,60} , {(ai,bi,ci) = (3, 1, 1), i = 61,...,75} .

Контейнер имеет объем F(u*) = 4825,16 . Выводы

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

Эффективность алгоритма подтверждается результатами для тестовых примеров.

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

Литература: 1. Uhler C., Wright S. J. Packing Ellipsoids with Overlap. SIAM Review, 55(4):671-706. 2013. 2. Kallrath J. Packing ellipsoids into volume-minimizing rectangular boxes. Journal of Global Optimization. D0I:10.1007/s10898-015-0348-6. 3. Pankratov A., Romanova T., Khlud O. Quasi-phi-functions in packing of ellipsoids. Radioelectronics & Informatics, 68:37-42. 2015.4. Lubachevsky B. D., Stillinger

F. H. Geometric properties of random disk packings. Journal of Statistical Physics, 60(5-6):561-583. 1990. 5. Bennell J.A., Oliveira J.F. A tutorial in irregular shape packing problem.Journal of the Operational Research Society. 2009. 60:93-105. 6. Chernov N., Stoyan Yu., Romanova T. Mathematical model and efficient algorithms for object packing problem // Computational Geometry: Theory and Applications. 2010. Vol. 43, № 5. P. 533-553. 7. Стецюк П.1., Романова Т.С., Субота 1.О. NLP-задача упаковки гомоте-тичних елшмв у прямокутний контейнер // Теорiя опти-мальних ршень. 2014. С. 139-146. 8. Stoyan Yu.G. A mathematical model and a solution method for the problem of placing various-sized circles into a strip / Yu.G. Stoyan,

G.N. Yaskov // European Journal of Operational Research. 2004. Vol. 156. P. 590-600. 9. Stoyan Y, Pankratov A, Romanova T. Quasi-phi-functions and optimal packing of ellipses. Journal of Global Optimization. 2015. DOI: 10.1007/ s10898-015-0331-2.

Поступила в редколлегию 23.09.2015

Рецензент: д-р техн. наук, ан.с. Гиль Н.И.

Хлуд Ольга Михайловна, аспирантка Института проблем машиностроения им. А.Н. Подгорного НАН Украины, Харьков. Научные интересы: исследование операций, математическое моделирование, геометрическое проектирование. Адрес: Украина, 61204, Харьков, пр-т Л.Свободы, д.36а, кв. 126, тел. +38(098)6166783.

Суббота Ирина Александровна, канд. техн. наук, преподаватель ХНУРЭ. Научные интересы: исследование операций, математическое моделирование, геометрическое проектирование. Адрес: Украина, 61183, Харьков, ул. Дружбы Народов, д.241, кв. 45, тел. +38(097) 250 36 17.

Романова Татьяна Евгеньевна, д-р техн. наук, профессор, ведущий научный сотрудник отдела математического моделирования и оптимального проектирования Института проблем машиностроения им. А.Н. Подгорного НАН Украины, Харьков. Научные интересы: исследование операций, математическое моделирование, геометрическое проектирование. Адрес: Украина, 61084, Харьков, ул. Новгородская, д.6а, кв. 31, тел. (057) 7013477.

Olga M. Khlud received Master's degree in System Analysis (2015) from Kharkiv National University of Radioelectronics. She is an aspirant at the Institute for Mechanical Engineering Problems of the National Academy of Sciences of Ukraine (Kharkiv). Her current research interests include mathematical modeling, operational research, packing and cutting. Address: Ukraine, 61204, Kharkiv, L.Svobody av, 36а, apt. 126, tel. +38(098)6166783

Irina A. Subota received Candidate of Technical Sciences degree in Mathematical Modeling and Computational Methods (2015) from Institute for Problems in Machinery of National Academy of Sciences of Ukraine (Kharkiv). From 2015 she is a teacher at the Department of System Engineering, Kharkiv National University of Radioelectronics. Her current research interests include mathematical modeling, operational research, packing and cutting.. Address: Ukraine, 61183, Kharkiv, str. Druzhby Narodov, 241, apt. 45, tel. +38(097) 250 36 17.

Tatiana E. Romanova received Doctor of Technical Sciences degree in Mathematical Modeling and Computational Methods (2003) from Institute of Cybernetics of the National Academy of Sciences of Ukraine (Kiev). From 2002 he is a senior researcher at the Department of Mathematical Modeling and Optimal Design, Institute for Mechanical Engineering Problems of the National Academy of Sciences of Ukraine. From 2005 she is a professor at the Department of Applied Mathematics, Kharkiv National University of Radioelectronics. Her current research interests include mathematical modeling, operational research, computational geometry, optimisation, packing, cutting and covering. Address: Ukraine, 61145, Kharkiv, Novgorodskaya str., 6a, apt. 31, tel.(057) 7013477.