Эффективность методов локального поиска в задаче минимизации энергии плоского кристалла Текст научной статьи по специальности «Математика»

| Прикладные проблемы оптимизации

УДК 519.853.4

Горчаков А.Ю., Посыпкин М.А.

Федеральный исследовательские! центр «Информатика и управление» РАН, г. Москва, Россия

ЭФФЕКТИВНОСТЬ МЕТОДОВ ЛОКАЛЬНОГО ПОИСКА В ЗАДАЧЕ МИНИМИЗАЦИИ

ЭНЕРГИИ ПЛОСКОГО КРИСТАЛЛА

Аннотация

Работа посвящена решению задачи, состоящей в нахождении минимума энергии совокупности атомов фрагмента плоской кристаллической решетки. Приводится математическая постановка задачи. Для ее решения применены два класса методов: квазиньютоновские методы семейства L-BFGS и методы координатного спуска, не требующие производных. Экспериментально показано, что в зависимости от постановки эффективны методы одного или другого класса.

Ключевые слова

Двумерные кристаллы; методы оптимизации; быстрое автоматическое дифференцирование; графеноподобные структуры; равновесное состояние.

Gorchakov A.Ju., Posypkin M.A.

Federal Research Center Computer Science and Control of the Russian Academy of Sciences, Moscow, Russia

THE EFFECTIVENESS OF LOCAL SEARCH METHODS IN THE PROBLEM OF FINDING THE

MINIMUM ENERGY OF A 2-D CRYSTAL

Abstract

The work is devoted to solving the problem of finding the energy minimum of the aggregate of atoms of a fragment of a planar crystal lattice. A mathematical statement of the problem is given. For its solution, two classes of methods are used: quasi-Newton methods of the L-BFGS family and coordinate descent methods that do not require derivatives. It has been shown experimentally that, depending on the formulation, methods of one or another class are effective.

Keywords

Two-dimensional crystals; optimization methods; fast automatic differentiation; graphene-like structures; equilibrium state.

1. Введение

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

соответствуют стабильным состояниям вещества.

Для вычисления энергии применяется метод,

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

дифференцируемость целевой функции.

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

2. Постановка оптимизационных задач

Рассматриваемая в даннои работе плоская модель многослойного кусочно-однородного материала предложена в [1, 2]. Материал представляется в виде периодической кусочно-однородной многослойной структуры, в рамках которои могут отличаться типы атомов в различных слоях. Модель накладывает ряд ограничении на структуру слоев:

• каждьш слои состоит из одинаковых атомов, в разных слоях могут находиться различные атомы;

• расстояния между соседними атомами в одном слое одинаковы, но в разных слоях они могут отличаться;

• в рассматриваемой системе выделяется

группа из К параллельных слоев, которая периодически повторяется в направлении

оси у ;

• число атомов в каждом слое и общее число слоев считается потенциально неограниченным.

На Рис. 1 представлен пример модели, в которои повторяется группа из трех слоев. Каждыи слой состоит из атомов своего вида.

.У

dr.

К

о о о о о о о о о о о о о о о о

<> 0 0 0 0 0 0 0.....

слое

i

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

Рассматривается две постановки (далее -постановка I и постановка II). В обоих случаях

минимизируется энергия

E (х)

совокупности

K

Рис. 1 Двумерная модель материала

В рамках даннои модели положение атомов определяется следующими параметрами:

И., 1 = 1,...,К г - расстояние между слоем с

номером 1 и предыдущим слоем;

d., 1 = 1,... ,К

г - смещение первого атома в слое

1 с положительной абсциссой относительно нулевои отметки;

e., 1=1,..., К

расстояние между атомами в

Совокупность значении перечисленных

атомов, расположенных на соседних слоях.

Параметрами оптимизационной задачи являются

x=(h1 ,d1 ,s1 ,...,hK,dK,sK) переменные l 1 1 1 K K K ,

составляющие вектор размерности n = 3 K . Таким образом, решаемая задача оптимизации относится к числу так называемых задач с параллелепипедными ограничениями и математически формулируется следующим образом:

E(хmin,l,<х,<u.,i=1,... ,n

i i i . (1)

В постановке I в моделируемую совокупность

входят атомы, абсцисса которых лежит в заданном

интервале [ 0 ] . При этом в зависимости от значении параметров, задающих конфигурацию, число и состав атомов, входящих в моделируемую совокупность, могут меняться. В постановке II считается, что атомы, входящие в эту совокупность, не меняются.

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

R 0

отсечения о , т.е. радиус, за пределами которого,

энергия взаимодействия атома с окружающими

атомами равна нулю. Это позволяет выразить

E (х) х=(х х )

энергию А атома с координатами 1' "2

конечнои суммои:

E(х) = Z v( х,У)

У eU (х)

г2:

где г - ■■ _ окрестность

радиуса Rо точки х , а v(х,у) - энергия взаимодействия атомов, находящихся в точках х и

у . В даннои статье рассматривается только случаи парного потенциала Леннарда-Джонса V(х,у)=||х-у||"12-2||х-у||-6 .

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

U (х ) = { у eR 2:||х-у||< RQ)

приводят к скачкообразному изменению энергии взаимодействия, вызванному выходом части атомов за границы фрагмента.

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

Л п

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

Сглаживающей назовем дважды

дифференцируемую функцию и( х) , такую что и(-1 )=1, и(1 ) = 0, и (-1 ) = и"(-1 )=и'(1 ) = и"(1)=0, и (х)<п при -1<х< 1

Пусть радиус отсечения равен Я и у (х) -потенциал парного взаимодействия. Тогда

определим сглаженный потенциал х)

следующим образом

V(х), при х<Я — d,

у( х )={ и ( — - —+11 V (х), при R-d <х < Я,

^ d )

0, при х > Я.

Несложно показать, что функция х) дважды

дифференцируема и обращается в ноль при х—Я . В качестве сглаживающей функции в данной работе взят следующий полином:

и( х)=-0.1875 х5+0.625 х3-0.9375 х+0.5 .

3. Описание методов, использованных для решения оптимизационных задач

Рассмотрим два метода поиска локального

экстремума. Первый - Бройдена-Флетчера-Гольдфарба-Шанно с ограниченным

использованием памяти (L-BFGS) [3]. Итерационная схема данного метода:

x+i=xk-akHkgk,Hk+i=F (sk>yk>Hk) > Sk Xk+1 xk,yk gk +1-gk, где ak - шаг метода, определяемый с помощью линейного поиска, gk - это градиент функций, вычисленный на k-ом шаге, а Hk - аппроксймацйя матрйцы обратной к матрйце вторых пройзводных функцйй, а F- это процедура ее обновленйя.

Как мы вйдйм, метод L-BFGS требует уменйя вычйслять градйент мйнймйзйруемой функцйй. Градйент вычйслялся с помощью методйкй быстрого автоматйческого дйфференцйрованйя [4], прй этом прйменялся пакет программного обеспеченйя Adept [5]. Данный пакет позволяет автоматйческй дйфференцйровать функцйй, напйсанные на языках С й С++. Он йспользует подход перегрузкй операторов (operator overloading), поэтому прй его йспользованйй требуется очень небольшая модйфйкацйя кода. В коде функцйй требуется все актйвные переменные объявить спецйальным тйпом adouble; на практйке просто пройзводйтся пойск й замена double на adouble; далее вызывается спецйальная функцйя для вычйсленйя градйента.

Второй метод является модйфйкацйей йзвестного метода коордйнатного спуска (ACD). Суть метода состойт в том, что на каждом шаге выполняется сдвйг в одном йз направленйй вдоль коордйнатных осей. Еслй в полученной точке значенйе целевой функцйй меньше, чем в йсходной, то далее она выбйрается в качестве йсходной. В протйвном случае сдвйг выполняется в протйвоположном направленйй. Еслй последовательные сдвйгй в направленйй коордйнатного вектора й в протйвоположном не далй уменьшенйя целевой функцйй, то велйчйна шага по этому направленйю уменьшается путем

умноженйя на коэффйцйент а <1 . Также вычйсляется аппроксймацйя градйента функцйй, в направленйй которого выполняется лйнейный пойск. Итерацйй повторяются до тех пор, пока значенйе максймального шага по всем направленйям не станет меньше, чем заданная

велйчйна

£

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

Эксперйменты проводйлйсь на персональном компьютере с процессором Intel Core i7-4770, 3.4 GHz, 8Gb оператйвной памятй.

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

0.5<h¡< 1 ^!=0; 0<d¡<4; 0<я,.<4 = 1,... ,4 . Координаты точек генерировались с использованием равномерного случайного распределения в заданных границах.

Сравнивались методы Монте-Карло с различными методами поиска локального минимума:

- L-BFGS с параметрами: количество итерации не более 100, количество итерации линейного поиска не более 3 (ШFGS3);

- L-BFGS с параметрами: количество итерации не более 100, количество итерации линейного поиска не более 20 (ШFGS20);

- ACD;

Результаты работы методов приведены в

таблице 1. Сравнение показывает, что в целом методы CD и ACD при решении оптимизационной задачи в постановке I оказываются эффективней L-BFGS. По мнению авторов, это связано с тем, что целевая функция не является непрерывной Покоординатные графики целевой функции в окрестности локального минимума со значением -409.33 можно увидеть на рис.2.

Темпы улучшения значения целевой функции для методов ACD и lBFGS20 приведены на рис.3. ACD достаточно быстро (за 150 секунд) находит минимум со значением -453.84 и далее пытается его улучшить. lBFGS20 находит минимум со значением -443.02, делая в 16 раз больше итерации Монте-Карло.

Таблица 1. Результаты решения оптимизационной задачи в постановке I

Метод Количество итераций Монте-Карло Значение найденного минимума

lBFGS3 9228 -409.33

lBFGS20 2049 -443.02

ACD 128 -453.84

-360

-370

-380

-390

400

410

420

\ \

ч \ ___

: 1 ■■■ ■. • *

Рис.2 Покоординатные срезы значений функции в окрестности минимума

о

-50 -100 -150 -200 -250 -300 -350

100

200

300

400

500

600

V

400 450 -500

.........*.....«

Рис.3 Сравнение скорости сходимости IBFGS20 (пунктирная линия, черные квадраты) и ACD (штриховая линия, серые

треугольники)

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

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

Информация о финансовой поддержке

Выполнено при финансовой поддержке проекта НШ-8860.2016.1 программы государственной поддержки ведущих научных школ Россиискои Федерации, РФФИ в рамках научных проектов № 17-07-00510, 17-07-00493, программ Президиума РАН I.33 «Математические модели и средства для исследования экономических и физических процессов с использованием

высокопроизводительных вычислении» и I.5 «Интеллектуальные информационные технологии и системы».

Литература

1. Лурье С.А., Посыпкин М.А., Соляев Ю.О. Метод идентификации масштабных параметров градиентнои теории упругости на основе численных экспериментов для плоских композитных структур // International Journal of Open Information Technologies, 2015. Т. 3. № 6. С. 1-6.

2. Ю.Г. Евтушенко, С.А.Лурье, М.А.Посыпкин, Ю.О.Соляев. Применение методов оптимизации для поиска равновесных состоянии двумерных кристаллов // Журнал вычислительной математики и математической физики. 2016. Т. 56. № 12. С. 50-59.

3. J. Nocedal, "Updating quasi-Newton matrices with limited storage/'Mathematics of Computation 35 (1980) 773-782.

4. Аида-Заде К.Р., Евтушенко Ю.Г. Быстрое автоматическое дифференцирование / / Математическое моделирование, 1989. Т. 1, С. 121-139.

5. Hogan R. J. Fast reverse-mode automatic differentiation using expression templates in C++ //ACM Transactions on Mathematical Software (TOMS). - 2014. - Т. 40. - №. 4. - С. 26.

6. Ю.Г. Евтушенко, М. А. Посыпкин. Применение метода неравномерных покрытии для глобальной оптимизации частично целочисленных нелинейных задач // Журнал вычислительной математики и математической физики, 2011, том 51, № 8, с. 1376-1389.

7. Посыпкин М.А. Решение задач глобальной оптимизации в среде распределенных вычислении / / Программные продукты и системы. № 1. 2010. С. 23-29.

References

1. Lur'e S.A., Posypkin M.A., Soljaev Ju.O. Metod identifikacii masshtabnyh parametrov gradientnoj teorii uprugosti na osnove chislennyh jeksperimentov dlja ploskih kompozitnyh struktur // International Journal of Open Information Technologies, 2015. T. 3. № 6. S. 1-6.

2. Ju.G. Evtushenko, S.A.Lur'e, M.A.Posypkin, Ju.O.Soljaev. Primenenie metodov optimizacii dlja poiska ravnovesnyh sostojanij dvumernyh kristallov / / Zhurnal vychislitel'noj matematiki i matematicheskoj fiziki. 2016. T. 56. № 12. S. 50-59.

3. J. Nocedal, "Updating quasi-Newton matrices with limited storage/'Mathematics of Computation 35 (1980) 773-782.

4. Ajda-Zade K.R., Evtushenko Ju.G. Bystroe avtomaticheskoe differencirovanie // Matematicheskoe modelirovanie, 1989. T. 1, S. 121139.

5. Hogan R. J. Fast reverse-mode automatic differentiation using expression templates in C++ //ACM Transactions on Mathematical Software (TOMS). - 2014. - T. 40. - №. 4. - S. 26.

6. Ju.G. Evtushenko, M. A. Posypkin. Primenenie metoda neravnomernyh pokrytij dlja global'noj optimizacii chastichno celochislennyh nelinejnyh zadach / / Zhurnal vychislitel'noj matematiki i matematicheskoj fiziki, 2011, tom 51, № 8, s. 1376-1389.

7. Posypkin M.A. Reshenie zadach global'noj optimizacii v srede raspredelennyh vychislenij / / Programmnye produkty i sistemy. № 1. 2010. S. 23-29.

Поступила: 25.07.2017

Сведения об авторах:

Горчаков Андрей Юрьевич, кандидат физико-математических наук, ведущий математик отдела прикладных проблем оптимизации Вычислительного центра

им. А.А. Дородницына, Федеральный исследовательский центр «Информатика и управление» РАН, andrgor12@ gmail. com ;

Во второй! серии экспериментов решалась оптимизационная задача в постановке II, в качестве начального приближения была взята точка, соответствующая минимальному значению энергии, полученной в первои серии экспериментов, со значением минимума -453.84. Сравнивались методы ACD и lBFGS20. Результаты приведены в таблице 2. Результат ожидаемый, в случае непрерывнои и дважды дифференцируемой функции L-BFGS работает существенно быстрее.

Таблица 2. Результаты решения оптимизационной

задачи в постановке II

Метод Время (сек) Количество вычислений функции (и градиента для L-BFGS) Значение найденного минимума

lBFGS20 1.06 127 -456.44

ACD 9.13 15541 -456.44

Заключение

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

Посыпкин Михаил Анатольевич, доктор фйзйко-математйческйх наук, доцент, заведующйй отделом Вычйслйтельного центра йм. А.А. Дороднйцына, Федеральный йсследовательскйй центр «Информатйка й управленйе» РАН, mposypkin@gmail.com.

Note on the authors:

Gorchakov Andrei Yu., Candidate of Physical and Mathematical Sciences, Leading mathematician, Dorodnicyn Computing Center, Federal Research Center «Computer Science and Control» of Russian Academy of Sciences, andrgor12@ gmail. com;

Posypkin Mikhail A., Doctor of Physical and Mathematical Sciences, head of the department, Dorodnicyn Computing Center, Federal Research Center «Computer Science and Control» of Russian Academy of Sciences, mposypkin@gmail.com.