CC BY
87
НАУЧНОЕ ИЗДАНИЕ МГТУ ИМ. Н. Э. БАУМАНА
НАУКА и ОБРАЗОВАНИЕ
Эл № ФС77 - 4 8211. Государственная регистрация №042 1200025. ISSN 1994-0408
электронный научно-технический журнал
Глобальная оптимизация на основе гибридизации алгоритмов поиска гармонии и роя частиц # 07, июль 2014
DOI: 10.7463/0714.0717653
профессор, д.ф.-м.н. Карпенко А. П.1а, Печенина Т. В.2, Буланов В. А.1
УДК 519.6
1Россия, МГТУ им. Н.Э. Баумана 2ООО АСКОН Санкт-Петербург apkarpenko@mail.ru
В работе представлена авторская гибридизация алгоритмов поиска гармонии и роя частиц. Алгоритм призван устранить или уменьшить влияние недостатков известного алгоритма Global-best Harmony Search (GHS), предложенного в 2008 г. в работе Omran M.G.H., Mahdavi M. Описана программная реализация алгоритма и результаты исследования его эффективности при решении тестовых задачах глобальной оптимизации многомерных многоэкстремальных функций Растригина и Шекеля, а также практически значимой задачи оптимизации ферменной конструкции. Результаты исследования показали эффективность принятых алгоритмических и программных решений. Так, в результате решения последней задачи были получены результаты, лучшие аналогичных результатов, полученных другими авторами.
Ключевые слова: Глобальная оптимизация, гибридизация, алгоритм поиска гармонии, алгоритм роя частиц
Введение
В 80-х годах прошлого века началась разработка класса стохастических поисковых алгоритмов глобальной оптимизации, называемых поведенческими, интеллектуальными, метаэвристическими, вдохновлёнными (инспирированными) природой, роевыми, много-агентными, популяционными и т.д. Используем последний термин. Популяционные алгоритмы предполагают одновременную обработку нескольких вариантов решения задачи оптимизации, в то время как в классических («траекторных») поисковых алгоритмах в области поиска эволюционирует только один кандидат на решение этой задачи.
Опыт решения сложных задач глобальной оптимизации с помощью популяционных алгоритмов показывает, что применение одного такого алгоритма далеко не всегда является эффективным. Поэтому в настоящее время большое внимание уделяют гибридизации популяционных алгоритмов глобальной оптимизации. Гибридные алгоритмы объединяют
различные алгоритмы либо одинаковые алгоритмы, но с различными значениями свободных параметров. При этом эффективность одного алгоритма может компенсировать слабость другого [1 - 3].
Целями данной работы являются разработка гибридного алгоритма глобальной оптимизации, основанного на алгоритмах поиска гармонии и роя частиц; программная реализация алгоритма; исследование его эффективности на ряде известных тестовых задач, а также задаче размернойоптимизации ферменной конструкции.
Алгоритм поиска гармонии (HarmonySearchalgorithm, HS) предложен в 2001 Гимом (Z. W. Geem) [4]. Алгоритм вдохновлен процессом поиска музыкантами гармонии в музыке. Ситуации идеальной гармонии звуков алгоритм HSставит в соответствие глобальный экстремум в задаче многомерной оптимизации, а процессу импровизации музыканта -процедуру поиска этого экстремума.
Основными особенностями алгоритма HS являются его простота и возможность использования для решения задач как непрерывной, так и дискретной оптимизации. В силу этого алгоритм привлекает значительное внимание исследователей и быстро развивается как в теоретическом, так и в практическом планах. Канонический алгоритм HS показал себя эффективным при решении значительного числа различных прикладных задач оптимизации. Однако алгоритм не свободен от недостатков, наиболее существенными из которых являются низкая эффективность локального поиска и сложность подбора значений свободных параметров алгоритма. Известно большое число работ, в которых предлагаются различные модификации канонического алгоритма HS, имеющие целью преодоление его указанных недостатков.
В работе [5] предложен вариант алгоритма HS с динамической настройкой свободных параметров (Improved HS, IHS).
Модификация алгоритма поиска гармонии, получившая название Global-best HS (GHS), представлена в работе [6]. Суть модификации заключается в добавлении в алгоритм IHS некоторых механизмов, используемых алгоритмом оптимизации роем частиц (Particle Swarm Optimization, PSO) [7]. Улучшенный вариант алгоритма GHS, названный авторами Improved Global Harmony Search (IGHS), предложен в работе [8].
Алгоритм HS с само адаптирующимися свободными параметрами рассмотрен в работе [9]. Ещё один вариант само адаптирующегося алгоритмаЯ^ получивший название (Novel Self-adaptive HS (NHS), представлен в работе [10].
Хорошо известно, что популяционные алгоритмы в своих канонических вариантах эффективны в исследовании пространства поиска и отыскании в нем перспективных областей, но медленно сходятся к локальным экстремумам. Поэтому предложено большое число гибридизаций алгоритма HS с алгоритмами локального поиска. Так в работе [11] представлен гибридный алгоритм, в котором в качестве локального используется алгоритм последовательного квадратичного программирования (Sequential Quadratic Programming). В работе [12] рассмотрена гибридизация алгоритма HS с алгоритмом имитации отжига (Simulated Annealing). Модификация алгоритма HS, основанная на объединении
этого алгоритма с алгоритмом дифференциальной эволюции (Differential Evolution, DE) представлена в работе [13].
В пункте 1 приводим постановку задачи. Пункт 2 посвящён базовым алгоритмам HS и PSO. В пункте 3 представлена схема предложенного гибридного алгоритма, названного
HSPSO. Пункт 4 содержит результаты вычислительного эксперимента с разработанным алгоритмическим и программным обеспечением. В пункте 5. представляем постановку задачи размерной оптимизации ферменной конструкции и приводим результаты её решения с помощью алгоритма HSPSO. В заключении формулируем основные результаты работы и перспективы её развития.
1. Постановка задачи
Рассматриваем детерминированную задачу глобальной безусловной минимизации
min F (X) = F (X *) = F *, (1)
X eR"
где F(X) - скалярная целевая функция, F(X ) = F - её искомое минимальное значение, X = (х^, Х2,...$" ) - " -мерный вектор варьируемых параметров, R" - " -мерное арифметическое пространство. Полагаем, что фитнесс-функция (p(X) также подлежит минимизации.
Параллелепипед допустимых значений вектора варьируемых параметров имеет вид D = {x | х-< х < х+, i e[l: "]}œ R", (2)
- +
где xi , xi - заданные константы.
Заметим, что в вычислительной практике задачу глобальной условной оптимизации обычно сводят к задаче безусловной оптимизации, например, методом штрафных функций (п. 5).
2. Базовые алгоритмы
Канонический алгоритм HS [4, 14]. Музыкантам ставим в соответствие индивидов Sj, а оркестру - популяцию S = {s, i e [l: N]}, где N— число индивидов популяции. Аккорду, который берет музыкант S- в данный момент времени, сопоставляем значение вектора варьируемых параметров X-. Гармонию звуков формализует глобальный минимум фитнесс-функции p(X) . Совокупность текущих координат Xi, i e[1: N] всех индивидов популяции образует (N х " )-матрицу памяти гармоний (Harmony Memory, HM). Число строк этой матрицы N носит название размер памяти гармоний и обозначается HMS.
Схема канонического алгоритма HS включает в себя следующие основные шаги.
1) Инициализируем алгоритм.
2) Формируем вектор гармонии.
3) Выполняем пошаговую настройку вектора гармонии.
4) Обновляем матрицу НМ.
5) Завершаем итераций, если условие окончания итераций выполнено, возвращаемся к шагу 2 в противном случае.
Инициализация. Начальные значения компонентов вектораX¡, полагаем равномерно распределёнными в гиперпараллелепипеде Б:
Хи = х- + и1(0;1) (*; -х-) I е[1: Щ], ] е[1: п]. (3)
Здесь и (0; 1) - вещественное случайное число, равномерно распределённое в интервале [0; 1]. Совокупность так построенных векторов Х1 образует исходную матрицу
памяти гармоний
HM
Х1,1 Х1,2 X2,1 X2,2
V XN ,1 XN,2
X
1,n
X
2, n
X
N,n
Формирование вектора гармонии X' выполняем по следующим правилам. С вероятностью ^ в качестве компоненты х\ ■ вектора Х\ используем соответствующую компоненту случайного вектора X. из текущей матрицы НМ:
xz, j
HMhj, /1 = U1(1: N), j e[1: n].
Здесь U (1: N) - целое случайное число, равномерно распределённое в интервале [1 : N ]. Данная операция моделирует воспроизведение музыкантом какого-либо аккорда из его памяти гармоний и называется операцией случайного выбора (random selection).
С вероятностью (1 — <£h) в качестве компоненты x' ■ берём величину, сгенерированную по формуле вида (3). Операция адекватна ситуации формирования абсолютно случайного аккорда из доступного музыканту диапазона.
Свободный параметр алгоритма имеет смысл вероятности
использования памя-
ти гармоний (Harmony Memory Considering Rate) и обозначается HMCR.
Пошаговая настройка вектора гармонии. Если компонента x'i ■ вектора X[ выбрана из HM, то выполняем следующие действия.
С вероятностью £ изменяем эту компоненту по формуле
X, j = X, j + <AN(0;1) (4)
или с вероятностью (1 — £ ) оставляем без изменений. Здесь W±gn- дискретная случайная
величина, с равной вероятностью принимающая значения ±1. Формула означает, что компонента x' , с равной вероятностью получает случайное положительное или отрицательное приращение величиной bw N (0; 1). Параметр £ определяет вероятность изменения шага (Pitch Adjusting Rate) и носит наименование PAR. Параметр представляет собой величину шага, которому придаётся смысл полосы пропускания инструмента (Pitch Bandwidth Rate). В каноническом алгоритме HS этот параметр обозначают BWR.
Обновление матрицы HM выполняем по следующей схеме. Вычисляем значение фитнесс-функции <(X.'), соответствующее сформированному вектору X[. Если
«(X'f ) < <р(X ) , то худший вектор Xi памяти гармоний заменяем вектором X\. Век-
тор X находим из условия
тах «Д. ) = <Xi ).
/е[1: N ]
Завершение итераций. Если условие окончания итераций выполнено, то в качестве
'ь
решения задачи используем лучший вектор памяти гармоний Х{ , доставляющий наи
меньшее значение функции приспособленности:
тт <р(Х; ) = у(Х ).
/е[1: N ] ' ь
Рекомендуют следующие значения основных свободных параметров алгоритма: ИМ8 < 50; ИМСЯ=0,7 ^ 0,95; РЛЯ= 0,1^ 0,5. Два последних параметра используют для регулирования соотношения диверсификационных и интенсификационных свойств алгоритма - меньшие значения этих параметров диверсифицируют поиск (за счёт, конечно, уменьшения скорости его сходимости).
Канонический алгоритм Р80 [15]. В каноническом алгоритме PSO координаты
Т
частицы si на итерации I > 0 определят векторX / = (Хц,X/ 2,...Хг п) , а на итерации
^ + 1) - вектор X'= (хг-1, х' 2,... X/ п)Т, где Т - символ транспонирования. Начальные
координаты частицы s, равны X , (0); /е[1: N].
Итерации в алгоритмеPSO выполняем по схеме
X=X/ + V,
V = Ь1У!-+ ип (0; ьс ) ®(Х* - X,) + ип (0; ЬБ ) ^(х** - Хг). (5)
Здесь использованы следующие обозначения: ® - символ прямого произведения
т
векторов; Ь7, Ьс, Ь8 - свободные параметры алгоритма; V = V(0 = (V 1, V 2,... V п) -
(nxl) -вектор приращения координат частицы; Vt = Vt (t — 1); X* - вектор координат частицы £г-, соответствующий её наилучшему значению фитнесс-функции за время поиска [0 : t], то есть
min (р(X* (т)) = ^X*);
ге[0: t]
X * - вектор координат соседней с данной частицы с наилучшим значением приспособленности, то есть
min^X*) = ^(X- ), (3 4)(8)
j^N/
где N - множество номеров частиц, являющихся «соседями» данной частицы st; Un (a; b) - (nx 1) -вектор вещественных случайных чисел, каждое из которых равномерно распределено в интервале [a; b]; a < b.
Первое слагаемое в формуле (5) выполняет функции памяти частицы о её перемещении на предыдущей итерации и называется инерционной компонентой. Второе слагаемое в
той же формуле называют когнитивной компонентой, поскольку оно формализует тен-
*
денцию частицы s. вернуться из текущего положения X ■ в положение X* . Третье слагаемое формулы (3) именуют социальной компонентой. Компонента отражает стремление частицы s переместиться в направлении своей наиболее успешной соседней частицы.
Модифицированный алгоритм GHS [6, 8]. Основная идея модификации заключается в изменении в алгоритме HS правила генерации нового созвучия таким образом, чтобы это созвучие с заданной вероятностью PAR могло подражать наилучшему созвучию в HM. Такая модификация добавляет социальную модель поведения частиц алгоритма Р80в канонический алгоритм HS. Псевдокод модифицированного алгоритма генерации нового созвучия имеет следующий вид:
for i e[1: HM], j e[1: n]do
iftf(0; 1)<HMCR then
X, j = x,, j, k = U[1: HMS ] iftf(0; 1)<PARthen x\ j = x*l, l = U[1: n]
else
x, j = x—+ U (0; 1) (x+ — x— )
Здесь x* i - l-ая компонента лучшего созвучия текущей памяти гармоний.
Алгоритм GHS сохраняет важные преимущества алгоритма HS - простота структуры и невысокая вычислительная сложность. В то же время, данная модификация не меняет процедуру обновления HM, используемую в каноническом алгоритме HS, суть которой состоит в подражании новых гармоний наилучшему решению в HM. В результате, в про-
цессе поиска из HM исчезают все решения, ухудшающие значение целевой функции. В совокупности с взятой из алгоритма PS0 социальностью, это позволяет ускорить процесс сходимости алгоритма, но повышает риск «застревания» алгоритма в локальном минимуме, поскольку разнообразие созвучий в HM уменьшается.
3. Гибридный алгоритм HSPSO
Гибридный алгоритм HSPSO призван устранить или уменьшить влияние указанных недостатков алгоритма GHS. Суть предложенной модификации алгоритма GHS заключается в следующем изменении процедуры обновления HM: если новое сгенерированное созвучие лучше, чем случайно выбранное из HM, то новое созвучие включается в HM, а выбранное исключается из неё. Заметим, что подобная стратегия используется в генетических алгоритмах, когда в популяции с некоторой вероятностьюсохраняют особей с невысокой приспособленностью, позволяя тем самымсохранить разнообразие популяции и диверсифицировать поиск.
о UQPSO
В качестве критерия останова алгоритм HS использует стагнацию вычислительного процесса: если лучшее решение, достигнутое алгоритмом, изменяется в течение заданного числа итераций NI на фиксированную величину, меньшую или равную
S > 0 , то алгоритм завершает своё выполнение.
Схема алгоритма имеет следующий вид.
Шаг 1. Инициализация алгоритма.
Шаг 2. Инициализация HM - заполнение HM случайно сгенерированными по формуле (1) созвучиями, число которых равно размеру гармонической памяти HMS.
Шаг 3. Импровизация нового созвучия. Новое созвучие создаётся по схеме канонического алгоритма HS и состоит из просмотра HM, регулировки тона и случайного выборы.
Шаг 4. Обновление HM. Если новый вектор решения X' лучше, чем случайно выбранный из HM, то вектор X включаем в HM, а выбранный вектор исключаем из неё.
Шаг 5. Проверка критерия останова. Прекращаем вычисления, если имеет место стагнация вычислительного процесса (или если достигнуто максимально допустимое число итераций). В противном случае переходим к шагу 3.
Так же, как алгоритм GHS [6, 8], алгоритм HSPSOиспользует динамические параметры PAR, bw:
PAR= PARt) = PARm + PARmax — PARmn t; (6)
NImax
С f z.min
ln b
bw = bw (t) = bmaxexp
w
max
bw
NI
max
(7)
Здесь PARmm, PARmax, bW™1, Ь"ех - заданные константы; NImax - заданное максимально допустимое число итераций.
4. Программная реализация и исследование эффективности
Программная реализация алгоритма HSPSO выполнена на языке программирования С++ в среде разработки MicrosoftVisualStudio 2008. Основные функции, используемые программой, перечислены ниже:
VoidInitialize_Harmony_Memory () - заполнение памяти гармоний HM случайно сгенерированными созвучиями;
VoidFindBestinHM () - поиск лучшего решения в HM; VoidFindWorstinHM () - поиск худшего решения в HM; VoidImprovisenewHarmony () - импровизацию нового созвучия; VoidUpdate_HM () - обновление HM;
VoidHarmonySearch () - основная функция, реализующая собственно алго-
ttciPSO
ритм HS .
Используем мультистарт с числом запусков алгоритма, равным 30. В качестве критериев эффективности алгоритма рассматриваем следующие величины: A - оценка вероятности попаданий в область точного решения с заданной погрешностью SF = 0,001 (точность алгоритма); MF - оценка математического ожидания найденного алгоритмом минимального значения целевой функции F(X); MI - оценка математического ожидания
числа итераций до начала процесса стагнации; ME - оценка математического ожидания числа испытаний (вычислений значений целевой функции); оценки среднегоквадратично-го отклонения трёх последних величин Ор, Ода, Осоответственно.
Вычислительные эксперименты выполнены на компьютере с процессором IntelCore2 DuoProcessor T7200 и оперативной памятью объёмом 2,00 ГБ под управлением операционной системой Windows7.
Выполнено исследование зависимости эффективности алгоритма HSPSOот значений его свободных параметров PAR и HMS. В качестве тестовых использованы два ярких представителя многоэкстремальных функций - функции Растригина и Шекеля [16]. Рассмотрены размерности вектора варьируемых параметров Хэтих функций, равные 8, 16, 32.
Поскольку в алгоритме HSPSO параметр PAR является динамическим (см. формулу (6)), варьировалась верхняя граница этого параметра PARax, принимая значения 0,4,
0,65, 0,9. Нижняя граница PAIRim принята постоянной и равной 0,01. Исследование выполнено для размеров памяти гармоний HMS=25, 50, 100. В соответствии с рекомендациями работы [6], значение параметра HMCR принято равным 0,95. Значения параметра
b менялись в процессе итераций от b"m = 0,001 до b"ax = 0,01 по формуле (7). Если не оговорено иное, использованы следующие значения остальных свободных параметров:
t
максимально допустимое число итераций N.1' = 10000; число итераций до начала
процесса стагнации N1^^ =1000; соответствующая точность £ = 10-6.
4.1. Функция Растригина
Решение отыскивается в области
П = {Х | - 5 < < 5, г е[1: и]}.
Минимальное значение функции Растригина в этой области достигается в точке
* л
X = 0 и равно нулю.
Результаты исследования представлены в таблицах 1 - 3 и на иллюстрирующих их рисунках 1 - 4.
Таблица 1 - Результаты исследования эффективности алгоритма при варьировании величин РАИ^ ИМБ\ функция Растригина; N1 = 10000; П = 8
НМ8 25 50 100
0,4 0,65 0,9 0,4 0,65 0,9 0,4 0,65 0,9
МР 0,1 0,07 0,26 0,1 0,23 0,17 0,07 0,13 0,23
аг 0,3 0,25 0,5 0,3 0,49 0,45 0,25 0,42 0,49
Р'(Х') 0 0 0 0 0 0 0 0 0
М1 8610 8278 6908 8792 7883 7786 9870 8691 8260
От 1290 1401 1657 1181 1457 1568 291 1411 1687
МЕ 8635 8303 6933 8842 7933 7836 9970 8791 8360
Оце 1290 1401 1657 1181 1457 1662 291 1411 1687
М\\Х-Х'\\ 0,08 0,07 0,28 0,11 0,17 0,15 0,07 0,11 0,17
А, % 90 90 76 90 80 86 90 90 80
Таблица 2 - Результаты исследования эффективности алгоритма при варьировании величин РАИ^, ИМБ\ функция Растригина; N1 = 10000; П = 16
ЯЛ/5 25 50 100
Р АРтпг 0,4 0,65 0,9 0,4 0,65 0,9 0,4 0,65 0,9
МР 0."2 0,62 0,75 0,69 0,79 1,1 2,58 1,10 1,17
аР 0,82 0,75 0,84 0,53 1,01 1,02 0,89 0,95 0,95
Р\Х') 0 0 0 0,01 0 0 1,08 0 0
М1 9996 9848 9701 10000 9964 9874 10000 10000 10000
От 20 572 643 0 193 471 0 0 0
МЕ 10021 9873 9726 10050 10014 9924 10100 10100 10100
Оце 20 572 643 0 193 471 0 0 0
М\\Х-Х'\ 1 0,47 0,45 0,4 0,7 0.58 0,74 0,86 0.66 0,87
А, % 36 43 40 0 26 20 0 0 20
Таблица 3 - Результаты исследования эффективности алгоритма при варьировании величин РЛЦ HMS: функция Растригина; N1 =10000; П = 32
HMS 25 50 100
p ARmax 0,4 0,65 0,9 0,4 0,65 0,9 0,4 0,65 0,9
MF 14,08 11,66 9,8 18,66 15,52 12,8 23,03 19,57 15,65
<TP 2,65 1,94 1,79 2,86 2,12 2,5 3,64 3,25 2,66
F'(X') 9,98 6,5 6,41 13,48 12,17 6,65 16,18 12,11 10,78
Ml 10000 10000 10000 10000 10000 10000 10000 10000 10000
ONI 0 0 0 0 0 0 0 0 0
ME 10025 10025 10025 10050 10050 10050 10100 10100 10100
0 0 0 0 0 0 0 0 0
M\\X-X-\\ 2,9 2,64 2,64 3,14 3,04 2,87 3,65 3,27 3,04
A, % 0 0 0 0 0 0 0 0 0
На указанных рисунках и далее, если не оговорено иное, пунктиром показаны величины МР , ММЕ ±&Е и т.д.
0,7 PAR 0,9
max
Рисунок 1 - Значения критерия эффективности MF в функции параметра PA: функция Растригина;
NT = 10000; HMS = 25; n = 8
max ' '
Рисунок 2 - Значения критерия эффективности МЕ в функции параметра ИМБ\ функция Растригина;
Л1тах = 10000; РАИах = 0,65; п = 8
Рисунок 3 - Значения критерия эффективности MF в функции параметра П : функция Растригина;
NImax = 10000; PARmax = 0,65 ;HMS=25
100% А 90% 80% 70% 60% 50% 40% 30% 20% 10% 0%
0 5 10 15 20 25 30 и 35
Рисунок 4 - Значения критерия эффективности A в функции параметра П : функция Растригина;
Ы1тах = 10000; РЛДпах = 0,65 -^=25
Из представленных результатов следует, что алгоритм с высокой вероятностью локализует глобальный минимум функции Растригина при размерности пространства поиска п=8. При этом наилучший результат по точности достигается при сочетании значений параметров HMS=25, РЛДпах = 0,65. С увеличением размерности пространства поиска, вероятность локализации минимума уменьшается, а среднее число итераций увеличивается до своего максимально допустимого значения N.I= 10000, то есть при высокой
размерности вектора X состояние стагнации вычислительного процесса не наступает.
Результаты вычислительных экспериментов, полученные при увеличении максимального числа итераций до N1 тах = 50000, представлены в таблице 4 и на рисунках 5,
6. Сплошная линия на рисунке 6 соответствует случаю, когда в качестве условия окончания итераций используется стагнация вычислений; пунктирная линия - достижение максимального числа итераций N1 тах.
Таблица 4 - Результаты исследования эффективности алгоритма при варьировании величины п: функция Растригина; М1тах = 50000; РЛ^Х = 0,65, ИМБ=25
п 8 16 32
0,07 0 0,01
мг 0,25 0,3 0.574
(Тр 0 0,637 0,8
8278 12474 27537
1401 2237 3862
8303 12499 27562
&КЕ 1401 2237 3862
М[\Х-Х*\) 0,07 0.243 0.467
Л, % 90 80 16
Рисунок 5 - Значения критерия эффективности МЕ в функции параметра П : функция Растригина;
= 50000; РЛЦ^х = 0,65 ИМ8=25
Рисунок 6 - Значения критерия эффективности A в функции параметра П : функция Растригина;
^тах = 50000; РЛ^ах = 0,65 ^=25
Результаты исследования показывают, что при увеличении максимально допустимого числа итераций точность алгоритма возрастает. Важно, что даже большое значение величины N1^^ = 1000 и малое значение величины £ = 10"6, не обеспечивают достижение состояния стагнации вычислительного процесса. Из рисунка 6 следует, что при использовании в качестве критерия окончания вычислений условия достижения максимально допустимого числа итераций NIтах, вероятность локализации глобального минимума функции Растригина кардинально повышается, так что даже при высокой размерности вектора варьируемых параметров алгоритм обеспечивает высокую надёжность.
4.2. Функция Шекеля
Поиск решения осуществляем в области
П = {х |0 < < 10, /е[1: п]}.
Число локальных минимумов функции Шекеля в этой области задано рав-
*
ным 10.Значение функции в точке глобального минимума равно ) = —10,1.
Результаты исследования представлены в таблицах 5 - 7 и на иллюстрирующих их рисунках7 - 10.
Таблица 5 - Результаты исследования эффективности алгоритма при варьировании величин РЛЕПах, НЖ: функция Шекеля; NIтах = 10000; П = 8
НМ8 25 50 100
РАЕтдЛ 0,4 0,65 0,9 0,4 0,65 0,9 0,4 0,65 0,9
№) -юл -10,1 -10,1 -10,1 -10,1 -10,1 -10,1 -10,1 -10,1
МР -4,94 -5,39 -5,7 -7,59 -6,5 -8,19 -8,54 -8,03 -8,05
3,86 4,14 4,44 2,77 4,05 3,58 2,75 3,01 3,52
М1 7060 6495 6205 9212 7866 7562 9934 9641 8782
#N1 1792 950 918 756 1178 1216 254 533 895
МЕ 7085 6520 6230 9262 7916 7612 10034 9741 8882
01/И 1792 950 918 756 1178 1216 254 533 895
м\\х-х*\\ 5,84 4,11 3,47 3,57 3,45 1,2 2,34 2,22 2
А, % 33 43 50 46 53 76 53 66 73
Таблица 6 - Результаты исследования эффективности алгоритма при варьировании величин НШ: функция Шекеля; N.I= 10000; П = 16
ЯМЗ 25 50 100
РА^тах 0,4 0,65 0,9 0,4 0,65 0,9 0,4 0,65 0,9
■ 10,1 -юд -10,1 -9,72 ■9.9 -10,1 -8.4 -9,1 ■ 10,1
МР -7,6 ■9,07 -8,05 -8,45 -9,56 -10,05 -7,42 -8,33 ■9,44
Ор 3,75 2,57 337 1,59 0,38 0,02 1 0,75 137
М1 10000 9990 9825 10000 10000 9996 10000 10000 10000
От 0 41 388 0 0 22 0 0 0
МЕ 10025 10015 9850 10050 10050 10046 10100 10100 10100
0 41 388 0 0 22 0 0 10000
М[ РГ-Х-П 1Р1 0,85 0,3 0,15 0.07 0.01 0.19 0,14 039
А, % 13 63 73 0 0 76 0 0 6
Таблица 7 - Результаты исследования эффективности алгоритма при варьировании величин НШ:функция Шекеля; = 10000; П = 32
НМБ 25 50 100
РАЕ 0,4 0,65 0,9 0.4 0,65 0,9 0,4 0,65 0,9
ЯГ) ■435 -4,57 -52 -232 -3,74 -3,94 -1,19 -1,72 -2 Д8
МР -2,81 -3,26 -4,12 -1,69 -236 -2,74 -0,97 -137 -1,85
аР 0,64 0,81 1,08 0,43 0,7 0,66 0.22 035 0,43
мл 10000 10000 10000 10000 10000 10000 9955 10000 10000
0 0 0 0 0 0 240 0 0
МЕ 10025 10025 10025 10050 10050 10050 10055 10100 10100
<?ПЕ 0 0 0 0 0 0 240 0 0
тх-х'\) 0,52 0,47 039 0,73 3,61 0,53 1 0,82 0,69
А, % 0 0 0 0 0 0 0 0 0
3 0 4 0 5 0 6 0, 7 0 8 0. 9 1Я та X
и
^ 1
и-" № ч
'-к' * * "1
л
г п
— *
1
Рисунок 7 - Значения критерия эффективности MF в функции параметра : функция Шекеля;
N1 х = 10000; ИЫБ = 50; П = 8
0 20 40 60 ЗОЯМ? 100
и
■ч
ч
(> — — - — <
\
ч
ч
ч
ч
•
1
-V.
Рисунок 8 - Значения критерия эффективности MF в функции параметра ИЫБ : функция Шекеля;
^тах = 10000; РА^х = 0,9 ; П = 8
10500
МЕ
10000 9500 9000 8500 8000 7500 7000 6500 6000
0 5 10 15 20 25 30 п 35
Рисунок 9 - Значения критерия эффективности МЕ в функции параметра П : функция Шекеля;
Мтах = 10000; РЛЦ^х = 0,9; НМБ = 50
100%
л
90% 80% 70% 60% 50% 40% 30% 20% 10% 0%
0 5 10 15 20 25 30 п 35
Рисунок 10 - Значения критерия эффективности А в функции параметра П : функция Шекеля;
Л^х = 10000; РЛ^ах = 0,9 ; НМ = 50
Из представленных результатов можно сделать выводы, аналогичные выводам предыдущего подраздела. Алгоритм НБР8° с высокой вероятностью локализует глобальный минимум функции Шекеля при размерности пространства поиска п=8. При этом наилучший результат по точности достигается при НМ5=50, □ = 0,9, и среднее число итераций алгоритма остаётся меньшим, чем заданное максимальное число итераций Шшах=10000. С увеличением размерности пространства поиска вероятность локализации
минимума уменьшается, а среднее число итераций увеличивается до своего максимального значения, то есть состояние стагнации вычислений не достигается.
Результаты экспериментов, соответствующие максимальному числу итераций Шшах=50000, представлены в таблице 8 и на соответствующих рисунках 11, 12. Сплошная линия на рисунке 12 соответствует случаю, когда в качестве условия окончания итераций используется стагнация вычислений; пунктирная линия - достижение максимального числа итераций ЖТтах.
Таблица 8 - Результаты исследования эффективности алгоритма при варьировании величины П : функция Шекеля; Мтах = 50000; РЛ^Х = 0,9; ИМ8=50
п 8 16 32
-10Д -10,065 -10,034
М¥ -3,19 -9,569 -10,031
3,58 1,359 0,001
М[№] 7562 14694 35295
Ш6 1516 2736
М [АТ] 7612 14744 35345
1216 1516 2736
1,2 0,004 0,006
А,% 76 56 30
40000 МЕ 35000
30000
25000
20000
15000
10000
5000
о -1———1—1————1————————————
0 10 20 30 п 40
Рисунок 11 - Значения критерия эффективности МЕ в функции параметра П : функция Шекеля; РЛ^Х = 0,9; ИМ8 = 50 ;М1тах=50000
Рисунок 12 - Значения критерия эффективности А в функции параметра П : функция Шекеля; РА^х = 0,9; НШ = 50 ; Мша*=50000
Как и в случае функции Растригина, состояние стагнации вычислительного процесса при высокой размерности вектора X не достигается. При использовании в качестве критерия окончания вычислений условия достижения максимально допустимого числа итераций Жтах, вероятность локализации глобального минимума рассматриваемой функции Шекеля повышается до 100%.
В целом, результаты исследования показывают, что для задач Растригина и Шекеля наилучшие значения свободных параметров различны. Так для функции Растригина наилучшие результаты были получены при ИМ£=25, = 0,65, а для функции Шекеля
- при ИМ£=50, РА^Х = 0,9. В обоих случаях результаты, полученные при использовании большого размера памяти гармоний ИМ£=100, не превзошли результатов, полученных сменьшими размерами этой памяти. Следовательно, небольшие размеры памяти гармоний предпочтительнее, поскольку обеспечивают экономию памяти ЭВМ. В алгоритме
применяется динамическая настройка параметра РАЯ(п. 3). При больших значениях верхней границы этого параметра, □ алгоритм показывает лучшие результаты.
При этом также уменьшается число итераций, необходимых для отыскания глобального минимума. Эффект объясняется тем, что в этом случае на ранних этапах поиска любое-созвучие, хранящееся в гармонической памяти, может участвовать всоздании нового решения, что обеспечивает диверсификацию поискового процесса. На более поздних итера-
циях новое созвучие сбольшей вероятностью подражает лучшему из созвучий памяти гармоний, тем самым ускоряя сходимость алгоритма.
5. Параметрическая оптимизация веса ферменной конструкции с помощью алгоритма
н8Р80
5.1. Постановка задачи
Рассматриваем задачу размерной оптимизации ферменной конструкции, состоящей из десяти балок, схема которой представлена на рисунке 13 [17].
Рисунок 13 - Ферменная конструкция из 10 балок
Целевая функция имеет смысл массы конструкции и определяется выражением
10
/ (х) = р£ л111,
I=1
где р- плотность материала конструкции; Д - площадь поперечного сечения 1-го элемента конструкции; / - длина этого элемента.
Ограничения накладываются на перемещения тх , Г узлов по осям 0 х 0у, а также на возникающие в балках напряжения (Гк и площади поперченного сечения Д. Перемещения в узлах 1, 5 должны равняться нулю: г = 0; г = 0; Г = 0; г = 0. Для остальных узлов имеют место ограничения типа неравенств вида
гдоп - rX] > 0, Гдоп -ryj > 0; j = 2,3,4,6,
где гои - заданное максимально допустимое смещение узла. Ограничения на напряжения в элементах конструкции и площади поперечных сечений балок имеют, соответственно, вид
(доп > 0,
Aax - A > о, A - Amin ^ о.
Здесь (Jdon- заданное максимально допустимое напряжение; 4т, Aтах- заданные допустимые площади поперечных сечений; i е [1:10].
Принимаем следующие соглашения: допустимая величина напряжений в стержнях (Гдоп = 172,434, Мпа; допустимые смещения узлов гои =0,0508, м; минимальная и макси-
2 2
мальная площади поперечных сечений балок равны 4т =0,000645, м , Aтах=0,03, м ; длины балок /=9,1438, м; силы P, P2равны P = 667,233, кН, P2 = 222,411, кН. Модуль Юнга материала конструкции принимаем равным Е=70,Гпа; плотность материала -
р = 268999, кг.
Погрешности выполнения ограничений для напряжений б и перемещений s принимаем равными Б( = 0,1, ( п=17,2434, Мпа и Бг = 0,1, Гоп= 0,005, м соответственно.
Нормальные напряжения в сечениях балок рассчитываем по формуле
N (i = ,
i 4 '
где N - нормальная сила, возникающая в i-ой балке под действием внешних сил; i е[1:10]. Полаем, что для удлинений Д^ балок справедлив закон Гука, так что
N /■
Д/{ = , i е[1:10]. E 4 L J
Перемещения каждого из узлов по осям 0x, 0y определяем как сумму удлинений
балок, соединяющихся в этом узле [17].
Для сведения поставленной задачи условной оптимизации к задаче безусловной оптимизации вида (1), (2), используем метод штрафных функций. В качестве F(X) используем функцию
K L
F(X) = f(X) + min(0, gk(X)) + Jdßl min(0, h(X)),
k=1 l=1
где (Хк, - штрафные коэффициенты; К, Ь - числа ограничений типа неравенств ^ (X) > 0 и равенств Ц (X) = 0 соответственно.
5.2. Вычислительный эксперимент
При решении поставленной задачи использованы следующие значения параметров алгоритма ШРЮ: #М5=30; ИМСЯ=0,95; РАРЩ^ = 0,1; □ РЛЦпал = 0,9; Ьтт = 0,0001; Ьшах = 0,001; N1Шах = 10000.
Уу ' 7 У ' -3 Шал
Полученное решение представлено в таблице 9. Соответствующие значения напряжений, смещений узлов и относительных погрешностей выполнения ограничений приведены в таблицах 10 - 14. Сходимость итерационного процесса иллюстрирует рисунок 14.
Таблица 9 - Результаты решения задачи
Величина поперечного сечения балки, м2. Масса, кг
Аг А2 А3 л4 Аб а7 ¿9 ¿10
0,013 0,0012 0,009 0,003 0,0036 0,0088 0,01 0,0037 0,006 0,001 2014,4
Таблица 10 - Напряжения в балках при найденных поперечных сечениях
Напряжения в балках,МПа.
0"! <?2 (Т4 "б а7 аэ ^10
73,2 -92,1 170 111 164 -161 -154 -150 -131 171
Относительные погрешности, %
58 47 1 35 5 7 11 13 24 1
Таблица 11 - Смещения узлов по оси 0 X при найденных сечениях
Смещения узлов
№ узла 1 2 3 4 5 6
гх, м. 0 0,0147 0,0490 -0,0006 0 0,0013
Погрешность, % 0 71 3,5 99 0 97
Таблица 12 - Смещения узлов по оси 0у при найденных сечениях
Смещения узлов
№ узла 1 2 3 4 5 6
гу,м. 0 -0,0150 -0,0053 0,0183 0 -0,0008
Погрешность, % 0 71 90 64 0 98
Рисунок 14-Сходимость итерационного процесса: Р - масса конструкции
Представленные результаты показывают, что в процессе оптимизации масса конструкции уменьшилась с 6500, кг до 2014,4 кг. Напряжения, возникающие в элементах конструкции при найденных площадях поперечных сечений балок, а также смещения узлов не превосходят своих допустимых значений.
Рассматриваемая задача оптимизации веса ферменной конструкции решалась различными методами [18]. Представленные результаты решения задачи показывают, что ал-
ттс<РЯО
горитм Но дает лучший результат по сравнению с результатами, полученными в предыдущих работах. Отметим, что решение, полученное алгоритмом , также превосходит решение, которое обеспечивает канонический алгоритм поиска гармонии [18].
Заключение
В работе представлена авторская гибридизация алгоритмов поиска гармонии и роя
частиц, названная HSPSO. Выполнена программная реализация алгоритма и исследование его эффективности на двух тестовых многомерных многоэкстремальных задачах глобальной оптимизации, а также на практически значимой задаче размерной оптимизации ферменной конструкции. Результаты исследования показали эффективность принятых алгоритмических и программных решений. Так, в результате решения последней задачи были получены результаты, лучшие аналогичных результатов, полученных другими авторами.
В развитие работы авторы планируют исследование эффективности алгоритма
HSPSO на более широком наборе тестовых задач, параллельную реализацию алгоритма.
Список литературы
1. Wang X. Hybrid nature-inspired computation method for optimization. Doct. Diss. Helsinki University of Technology, TKK Dissertations, 2009. 161 p.
2. El-Abd, Kamel M. A taxonomy of cooperative search algorithms // In: Hybrid Metaheuristics. Proc. Second International Workshop. Springer Berlin Heidelberg, 2005. P. 32-41. (Ser. Lecture Notes in Computer Science; vol. 3636). DOI: 10.1007/11546245 4
3. Raidl G.R. A Unified View on Hybrid Metaheuristics // In: Hybrid Metaheuristics. Proc. Third International Workshop. Springer Berlin Heidelberg, 2006. P. 1-12. (Ser. Lecture Notes in Computer Science; vol. 4030). DOI: 10.1007/11890584 1
4. Geem Z.W., Kim J.H., Loganathan G.V. A new heuristic optimization algorithm: harmony search // Simulation. 2001. Vol. 76, no. 2. P. 60-68.
5. Mahdavi M., Fesanghary M., Damangir E. An improved harmony search algorithm for solving optimization problems // Applied Mathematics and Computation. 2007. Vol. 188, no. 2. P. 1567-1579. DOI: 10.1016/j.amc.2006.11.033
6. Omran M.G.H., Mahdavi M. Global-best harmony search // Applied Mathematics and Computation. 2008. Vol. 198, no. 2. P. 643-656. DOI: 10.1016/j.amc.2007.09.004
7. Карпенко А.П., Селиверстов Е.Ю. Обзор методов роя частиц (PSO) для задачи глобальной оптимизации» // Наука и образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2009. № 3. Режим доступа: http://technomag.edu.ru/doc/116072.html (дата обращения 01.06.2014).
8. Lin Y., Chang-Ming X. A new improved harmony search algorithm for continuous optimization problem // 2011 International Conference on Computer Science and Network Technology (ICCSNT). Vol. 2. IEEE, 2011. P. 686-690. DOI: 10.1109/ICCSNT.2011.6182059
9. Wang Chia-Ming, Huang Yin-Fu. Self-adaptive harmony search algorithm for optimization // Expert Systems with Applications. 2010. Vol. 37, no. 4. P. 2826-2837. DOI: 10.1016/j.eswa.2009.09.008
10. Chen Jing, Man Hong-Fang, Wang Ya-Min. Novel Self-adaptive Harmony Search Algorithm for continuous optimization problems // 2011 30th Chinese Control Conference (CCC). IEEE, 2011. P. 5452-5456.
11. Fesanghary M., Mahdavi M., Minary-Jolandan M., Alizadeh Y. Hybridizing harmony search algorithm with sequential quadratic programming for engineering optimization problems// Computer Methods in Applied Mechanics and Engineering. 2008. Vol. 197, no. 33-40. P. 3080-3091. DOI: 10.1016/j.cma.2008.02.006
12. Taherinejad N. Highly reliable harmony search algorithm // ECCTD 2009. European Conference Circuit Theory and Design. IEEE, 2009. P. 818-822. DOI: 10.1109/ECCTD.2009.5275109
13. Gao X.Z., Wang X., Ovaska S.J., He Xu. A modified harmony search method in constrained optimization // International Journal of Innovative Computing, Information and Control. 2009. Vol. 6, no. 9. P. 4235-4247.
14. Xin-She Yang. Harmony Search as a Metaheuristic Algorithm // In: Music-Inspired Harmony Search Algorithm: Theory and Applications / Z.W. Geem, ed. Springer Berlin Heidelberg, 2009, pp. 1-14. (Ser. Studies in Computational Intelligence; vol. 191). DOI: 10.1007/978-3-642-00185-7 1
15. Kennedy J., Eberhart R. Particle swarm optimization // Proceedings of IEEE International Conference on Neural Networks. Vol. 4. IEEE, 1995. P. 1942-1948. DOI: 10.1109/ICNN.1995.488968
16. Tang K., Yao X., Suganthan P.N., MacNish C., Chen Y.P., Chen C.M., Yang Z. Benchmark Functions for the CEC'2008 Special Session and Competition on Large Scale Global Optimization. Technical Report. Nature Inspired Computation and Applications Laboratory, USTC, China, 2007. Available at: http://nical.ustc.edu.cn/cec08ss.php , accessed 01.06.2014.
17. Hultman M. Weight Optimization of Steel Trusses by a Genetic Algorithm-Size, Shape and Topology Optimization According to Eurocode. Report TVBK. Department of Structural Engineering Lund Institute of Technology, Lund University, 2010. 61 p.
18. Lee K.S., Geem Z.W. A new meta-heuristic algorithm for continuous engineering optimization: harmony search theory and practice // Computer Methods in Applied Mechanics and Engineering. 2005. Vol. 194, iss. 36-38. P. 3902-3933. DOI: 10.1016/j.cma.2004.09.007
scientific periodical of thïï bauman mstu
SCIENCE and EDUCATION
EL № FS77 - 48211. N»0421200025. ISSN 1994-0408
electronic scientific and technical journal
Global Optimization Based on the Hybridization of Harmony
Search and Particle Swarm Optimization Methods
# 07, July 2014
DOI: 10.7463/0714.0717653
A.P. Karpenko1a, T.V. Pechenina2, V.A. Bulanov1
1Bauman Moscow State Technical University, 105005, Moscow, Russian Federation
2ASCON 199155, St-Petersburg, Russian Federation
apkarpenko@mail.ru
Keywords: global optimization, PSO, hybridization
We consider a class of stochastic search algorithms of global optimization which in various publications are called behavioural, intellectual, metaheuristic, inspired by the nature, swarm, multi-agent, population, etc. We use the last term.
Experience in using the population algorithms to solve challenges of global optimization shows that application of one such algorithm may not always effective. Therefore now great attention is paid to hybridization of population algorithms of global optimization. Hybrid algorithms unite various algorithms or identical algorithms, but with various values of free parameters. Thus efficiency of one algorithm can compensate weakness of another.
The purposes of the work are development of hybrid algorithm of global optimization based on known algorithms of harmony search (HS) and swarm of particles (PSO), software implementation of algorithm, study of its efficiency using a number of known benchmark problems, and a problem of dimensional optimization of truss structure.
We set a problem of global optimization, consider basic algorithms of HS and PSO, give a
ttoPSO
flow chart of the offered hybrid algorithm called , present results of computing experi-
ments with developed algorithm and software, formulate main results of work and prospects of its development.
References
1. Wang X. Hybrid nature-inspired computation method for optimization. Doctoral Dissertation. Helsinki University of Technology, TKK Dissertations, 2009, 161 p.
2. El-Abd, Kamel M. A taxonomy of cooperative search algorithms. In: Hybrid Metaheuristics. Proc. Second International Workshop. Springer Berlin Heidelberg, 2005, pp. 32-41. (Ser. Lecture Notes in Computer Science; vol. 3636). DOI: 10.1007/11546245 4
3. Raidl G.R. A Unified View on Hybrid Metaheuristics. In: Hybrid Metaheuristics. Proc. Third International Workshop. Springer Berlin Heidelberg, 2006, pp. 1-12. (Ser. Lecture Notes in Computer Science; vol. 4030). DOI: 10.1007/11890584 1
4. Geem Z.W., Kim J.H., Loganathan G.V. A new heuristic optimization algorithm: harmony search. Simulation, 2001, vol.76, no. 2, pp. 60-68.
5. Mahdavi M., Fesanghary M., Damangir E. An improved harmony search algorithm for solving optimization problems. Applied Mathematics and Computation, 2007, vol. 188, iss. 2, pp. 1567-1579. DOI: 10.1016/j.amc.2006.11.033
6. Omran M.G.H., Mahdavi M. Global-best harmony search. Applied Mathematics and Computation,, 2008, vol. 198, no. 2, pp. 643-656. DOI: 10.1016/j.amc.2007.09.004
7. Karpenko A.P., Seliverstov E.Yu. Review of the particle swarm optimization method (PSO) for a global optimization problem. Nauka i obrazovanie MGTU im. N.E. Baumana = Science and Education of the Bauman MSTU, 2009, no. 3. Available at: http://technomag.edu.ru/doc/116072.html , accessed 01.06.2014. (in Russian).
8. Lin Y., Chang-Ming X. A new improved harmony search algorithm for continuous optimization problem. 2011 International Conference on Computer Science and Network Technology (ICCSNT). Vol. 2. IEEE, 2011, pp. 686-690. DOI: 10.1109/ICCSNT.2011.6182059
9. Wang Chia-Ming, Huang Yin-Fu. Self-adaptive harmony search algorithm for optimization. Expert Systems with Applications, 2010, vol. 37, no. 4, pp. 2826-2837. DOI: 10.1016/j.eswa.2009.09.008
10. Chen Jing, Man Hong-Fang, Wang Ya-Min. Novel Self-adaptive Harmony Search Algorithm for continuous optimization problems. 2011 30th Chinese Control Conference (CCC). IEEE, 2011, pp. 5452-5456.
11. Fesanghary M., Mahdavi M., Minary-Jolandan M., Alizadeh Y. Hybridizing harmony search algorithm with sequential quadratic programming for engineering optimization problems. Computer Methods in Applied Mechanics and Engineering, 2008, vol. 197, iss. 33-40, pp. 3080-3091. DOI: 10.1016/j.cma.2008.02.006
12. Taherinejad N. Highly reliable harmony search algorithm. ECCTD 2009. European Conference Circuit Theory and Design. IEEE, 2009, pp. 818-822. DOI: 10.1109/ECCTD.2009.5275109
13. Gao X.Z., Wang X., Ovaska S.J., He Xu. A modified harmony search method in constrained optimization. International Journal of Innovative Computing, Information and Control, 2009, vol. 6, no.9, pp. 4235-4247.
14. Xin-She Yang. Harmony Search as a Metaheuristic Algorithm. In: Geem Z.W., ed. Music-Inspired Harmony Search Algorithm: Theory and Applications. Springer Berlin Heidelberg, 2009, pp. 1-14. (Ser. Studies in Computational Intelligence; vol. 191). DOI: 10.1007/978-3642-00185-7 1
15. Kennedy J., Eberhart R. Particle swarm optimization. Proceedings of IEEE International Conference on Neural Networks. Vol. 4. IEEE, 1995, pp. 1942-1948. DOI: 10.1109/ICNN.1995.488968
16. Tang K., Yao X., Suganthan P.N., MacNish C., Chen Y.P., Chen C.M., Yang Z. Benchmark Functions for the CEC'2008 Special Session and Competition on Large Scale Global Optimization. Technical Report. Nature Inspired Computation and Applications Laboratory, USTC, China, 2007. Available at: http://nical.ustc.edu.cn/cec08ss.php, accessed 01.06.2014.
17. Hultman M. Weight Optimization of Steel Trusses by a Genetic Algorithm-Size, Shape and Topology Optimization According to Eurocode. Report TVBK, Department of Structural Engineering Lund Institute of Technology, Lund University, 2010. 61 p.
18. Lee K.S., Geem Z.W. A new meta-heuristic algorithm for continuous engineering optimization: harmony search theory and practice. Computer Methods in Applied Mechanics and Engineering,, 2005, vol. 194, no. 36-38, pp. 3902-3933. DOI: 10.1016/j.cma.2004.09.007
