Журнал вычислительной математики и математической физики, 2023, T. 63, № 1, стр. 61-73

ВВЕДЕНИЕ

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

Можно выделить несколько различных направлений в разработке численных методов решения задач оптимального управления, которые существенно отличаются друг от друга. Прежде всего следует упомянуть прямые методы, основанные на спуске в пространстве управлений. Ряд исследований связан с непрямыми методами, в которых с помощью принципа максимума Л.С. Понтрягина исходная задача оптимального управления сводится к краевой (см. [1]). В [2] Р.П. Федоренко предложил подход, который базируется на использовании идей метода линеаризации. Большой цикл работ по численным методам оптимального управления был выполнен Н.Н. Моисеевым (см. [3]) и его учениками.

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

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

Второй способ вычисления градиента функционала – использование сопряженной задачи. Здесь представляется естественным сформулировать исходную задачу оптимального управления в “непрерывном” пространстве, сформулировать сопряженную задачу также в “непрерывном” пространстве, определить градиент функционала в “непрерывном” случае, и только после этого осуществлять аппроксимацию функционала, прямой и сопряженной задач. В результате этих действий получаются формулы для определения градиента целевого функционала в дискретном варианте. Однако здесь необходимо учитывать то, что получаемый дискретный градиент функционала неточен. Его значение зависит от способа дискретизации как прямой и сопряженной задач, так и от способа дискретизации самого градиента. Эта неточность может быть уменьшена, если уменьшить погрешности указанных аппроксимаций, или даже полностью ликвидирована за счет выбора “согласованных” аппроксимаций прямой и сопряженной задач и “согласованной” с ними аппроксимации градиента. Эти “согласованные” аппроксимации можно случайно угадать, но указать способ определения “согласованных” аппроксимаций, используя данный подход, по-видимому, невозможно. Если выбрать неправильный дискретный вариант сопряженной задачи, получается неверное значение градиента. С такими случаями приходилось встречаться на практике (см., например, [5]).

Другой путь связан с формулировкой рассматриваемой задачи оптимального управления в конечномерном пространстве и с определением градиента функционала сразу в дискретном виде. Если дискретный градиент определять с помощью метода множителей Лагранжа (сопряженных переменных), то такой путь известен как обобщенная методология быстрого автоматического дифференцирования (БАД-методология).

БАД-методология значительно превосходит по эффективности аналитическое дифференцирование (точнее, использование аналитически полученных формул) и вычисление производных с помощью метода конечных разностей. В работах Ю.Г. Евтушенко [6–8] предложен общий подход к дифференцированию сложных функций. Там показано, что БАД-методология дает возможность с единой позиции рассмотреть самые разнообразные задачи. Например, из общих формул дифференцирования, приведенных в этих работах, несложно получить формулы БАД-методологии для определения градиента функций многих переменных. При этом вычисление значения функции представляется как многошаговый процесс, в котором введены новые фазовые переменные. Эти фазовые переменные являются функциями независимых переменных, по которым ищутся производные заданной функции. Показано, что при определенных условиях использование БАД-методологии позволяет получить точное значение градиента дискретной задачи оптимального управления (дискретного целевого функционала при наличии дискретных связей) за время, не превосходящее утроенного времени вычисления самой функции.

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

В настоящей работе проводится сравнение разных подходов к вычислению градиента сложной функции многих переменных. Сравниваемыми подходами являются: использование точных, аналитически выведенных формул; использование формул, полученных с помощью БАД-методологии; использование стандартных программных пакетов, реализующих идеи БАД-методологии. Сравнение подходов осуществляется на примере сложной функции, представляющей энергию системы атомов, потенциал взаимодействия которых – потенциал Терсоффа.

1. СРАВНЕНИЕ РАЗЛИЧНЫХ МЕТОДОВ ДИФФЕРЕНЦИРОВАНИЯ ФУНКЦИЙ

Функция $f(u)$ называется элементарной функцией, если она может быть представлена как конечная суперпозиция основных элементарных функций и арифметических действий. Будем считать, что дифференцируемая скалярная функция $f(u)$ векторного аргумента $u \in {{R}^{r}}$ задана явно. Задача состоит в вычислении частных производных функции $f(u)$ по компонентам вектора $u$. В вычислительной практике часто встречаются функции, для которых точное вычисление частных производных представляет собой непростую задачу.

В таком случае можно воспользоваться несколькими способами нахождения частных производных функции $f(u)$_:

1. Приближенное численное дифференцирование. Если $f(u)$ – функция $r$-мерного вектора $u$, то при использовании самой простой формулы численного дифференцирования необходимо по меньшей мере $r + 1$ раз вычислить значение самой функции $f(u)$_, и $r$ раз произвести вычитание и деление. Результаты таких расчетов будут точными только для линейных функций; во всех остальных случаях градиент определяется с некоторой погрешностью. Ниже будет показано, что при вычислении частных производных энергии группы атомов (потенциал взаимодействия атомов – потенциал Терсоффа) по специфичным параметрам потенциала Терсоффа с помощью этого метода необходимо проводить дополнительные исследования, связанные с выбором приращения каждого параметра потенциала, причем при каждом его новом значении.

2. Аналитическое дифференцирование. Здесь имеется в виду использование выведенных аналитических выражений для вычисления частных производных функции $f(u)$. Следует также отметить, что в современной компьютерной литературе такой подход нередко называют прямым методом автоматического дифференцирования, или аналитическим дифференцированием. Если вектор $u$ – вектор большой размерности, то при применении этой техники возникают большие трудности.

3. Использование формул БАД (БАД-методология) (см. [6]–[8]). При использовании этого метода вычисление значения функции $f(u)$ записывается в виде многошагового процесса в каноническом виде, формулируется сопряженная задача в каноническом виде и с их помощью выводятся аналитические выражения для вычисления частных производных рассматриваемой функции.

4. Использование стандартных пакетов программ автоматического дифференцирования.

а) Используются пакеты программ для генерации кода вычисления производных по коду вычисления функции (source-to-source AD tool). В качестве примера такой программы можно привести TAPENADE (см. [9]).

б) Используются пакеты программ, перегружающих операторы в коде вычисления функции. В результате такой перезагрузки код вычисляет не только значение функции, но и значение ее производных (Operator Overloading AD tool). В качестве примера такой программы можно привести CoDiPack и Adept (см. [10], [11]).

Ниже в этом разделе приводится сравнение дифференцирования сложной функции с помощью формул БАД-методологии (способ 3) и с помощью аналитически полученных выражений для ее частных производных (способ 2).

Для того чтобы воспользоваться техникой БАД (способ 3), необходимо представить вычисление элементарной функции $f(u)$ как многошаговый процесс. С этой целью введем вектор $z \in {{R}^{k}}$ фазовых переменных (вектор вспомогательных, промежуточных значений функции):

где ${{Z}_{i}}$ – набор элементов ${{z}_{j}}$, которые входят в правую часть равенства ${{z}_{i}} = F(i,{{Z}_{i}},{{U}_{i}})$, а ${{U}_{i}}$ – набор элементов ${{u}_{j}}$, которые входят в правую часть этого равенства. Последовательность (1.1) строится таким образом, чтобы процесс был явным, и последняя компонента ${{z}_{k}}$ совпала бы со значением функции $f(u)$, т.е. ${{z}_{k}} = f(u)$.

Предполагаем, что все функции $F(i,{{Z}_{i}},{{U}_{i}})$ дифференцируемы по всем компонентам векторов $z$ и $u$. Введем вектор сопряженных переменных $p \in {{R}^{k}}$, который определим как решение следующей системы линейных алгебраических уравнений:

Тогда, согласно формулам БАД-методологии (см. [8]), компоненты градиента функции многих переменных $f(u)$ вычисляются по формуле

Таким образом, мы сначала последовательно выполняем вычисления в соответствии с $k$-шаговым процессом и запоминаем значения переменных $z$ и $u$ (прямой ход), а потом идем в обратную сторону и рассчитываем сопряженные переменные $p$ (обратный ход).

Что касается вычисления градиента функции методом 2, то здесь сопряженные переменные не определяются и не используются. Учитывая тот факт, что сложная функция $f(u)$ представляет собой суперпозицию более простых функций, представим вычисление ее значения в виде $k$-шагового алгоритма (1.1).

Тогда производные $\frac{{\partial f}}{{\partial {{u}_{m}}}}$, $m = \overline {1,r} $, вычисляются аналитически последовательно следующим образом:

Работу с формулами (1.2), (1.3) и (1.4) мы вначале проиллюстрируем на простом примере. Пусть задана элементарная функция $f(u) = f({{u}_{1}},{{u}_{2}}) = {{e}^{{{{u}_{1}} + {{{({{u}_{2}})}}^{2}}}}} + \sin {{\left( {{{{({{u}_{1}})}}^{3}} + {{u}_{2}}} \right)}^{2}}$. Для применения формул БАД-методологии (метод 3) строим многошаговый процесс ее вычисления следующим образом:

Используя формулы (1.2) и (1.3), имеем Сопряженные переменные вычисляются в обратном порядке: ${{p}_{6}}$, ${{p}_{5}}$, …, ${{p}_{1}}$.

Для того чтобы получить частные производные этой функции без использования сопряженной задачи (с помощью формул (1.4), метод 2), введем обозначения ${{\tilde {z}}_{i}} = \frac{{\partial {{z}_{i}}}}{{\partial {{u}_{1}}}}$, ${{\tilde {\tilde {z}}}_{i}} = \frac{{\partial {{z}_{i}}}}{{\partial {{u}_{2}}}}$. Значения ${{\tilde {z}}_{i}}$ и ${{\tilde {\tilde {z}}}_{i}}$ вычисляются по формулам

Тогда $\frac{{\partial f}}{{\partial {{u}_{1}}}} = {{\tilde {z}}_{6}}$, $\frac{{\partial f}}{{\partial {{u}_{2}}}} = {{\tilde {\tilde {z}}}_{6}}$.

В случае метода 2 вместо одного многошагового процесса (1.6), который требуется для вычисления сопряженных переменных, необходимо выполнить вычисления по $r = 2$ (количество независимых переменных) многошаговым процессам: (1.8) и (1.9). Следовательно, время, необходимое для вычисления градиента с помощью формул БАД, примерно в $r$ раз меньше времени, которое необходимо для вычисления этого градиента с помощью формул (1.4). Эффективность использования техники БАД возрастает, если растет количество независимых переменных.

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

где Здесь $I$ – количество атомов в рассматриваемой группе; $R$ и ${{R}_{{{\text{cut}}}}}$ – параметры, определяющиеся из полученных экспериментально геометрических характеристик вещества; ${{x}_{i}} = ({{x}_{{1i}}},{\kern 1pt} {{x}_{{2i}}},{\kern 1pt} {{x}_{{3i}}})$ – координаты $i$-го атома; ${{r}_{{ij}}} = \sqrt {{{{\left( {{{x}_{{1i}}} - {{x}_{{1j}}}} \right)}}^{2}} + {{{\left( {{{x}_{{2i}}} - {{x}_{{2j}}}} \right)}}^{2}} + {{{\left( {{{x}_{{3i}}} - {{x}_{{3j}}}} \right)}}^{2}}} $ – расстояние между атомами с номерами $i$ и $j$; ${{\theta }_{{ijk}}}$ – угол, образованный векторами, соединяющими атом $i$ с атомами $j$ и $k$ соответственно. Косинус этого угла вычисляется по формуле ${{q}_{{ijk}}} = \cos ({{\theta }_{{ijk}}}) = \frac{{r_{{ij}}^{2} + r_{{ik}}^{2} - r_{{jk}}^{2}}}{{2r_{{ij}}^{{}}r_{{ik}}^{{}}}}$.

Потенциал Терсоффа зависит от десяти параметров ($m = 10$), специфичных для моделируемого вещества (${{D}_{e}}$, ${{r}_{e}}$, $\beta $, $S$, $\eta $, $\gamma $, $\lambda $, $c$, $d$, $h$), и от координат атомов $x = {{({{x}_{1}},{{x}_{2}},...,{{x}_{I}})}^{T}}$. Указанные специфические параметры вещества, как правило, неизвестны. Для решения задач параметрической идентификации потенциалов межатомного взаимодействия в последнее время все большую актуальность приобретает применение различных методов оптимизации. В связи с этим возникает необходимость определения точного значения градиента функции энергии $E$ по специфическим для моделируемого вещества параметрам потенциала Терсоффа. Формулы для вычисления компонент этого градиента были получены в [12] с помощью БАД-методологии.

Оптимизационные задачи решаются при определенном, фиксированном положении базисных атомов рассматриваемой кристаллической структуры. Решив задачу идентификации параметров, нет уверенности, что положения базисных атомов будут соответствовать минимуму потенциальной энергии системы. Поэтому следующим шагом исследований является оптимизация по координатам частиц, расставляющая частицы в положения, соответствующие минимуму суммарной потенциальной энергии рассматриваемой системы атомов. На этом этапе возникает необходимость вычисления градиента функции энергии по координатам атомов. Формулы для вычисления этого градиента также были получены в [12] на основе формул БАД-методологии (метод 3).

Ниже мы построим вычислительный процесс для определения указанного выше градиента с помощью метода 2 (аналитического дифференцирования) и сравним его с тем, который получен с помощью метода 3.

Функция энергии $E$ представляет собой достаточно сложную суперпозицию более простых функций. Вычисление самой функции было представлено в [12] в виде некоторого многошагового процесса.

Пусть $\bar {u}$ и $\bar {z}$ – векторы, имеющие координаты

где

Каждая компонента ${{z}_{l}}$ вектора $\bar {z}$ состоит из множества других компонент, а процесс вычисления функции энергии $E$ строится следующим образом. Для всех фиксированных пар индексов $(i,j)$, таких что $i = \overline {1,I} $_, $j = \overline {1,I} $ и $j \ne i$, вычисляются значения $z_{{13}}^{{ij}}$. Далее для всех таких пар индексов $(i,j)$ и для всех индексов $k = \overline {1,I} $, таких что $k \ne i$ и $k \ne j$, вычисляются значения $z_{1}^{{ijk}},z_{2}^{{ijk}},...,z_{8}^{{ijk}}$. Компоненты $z_{1}^{{ijk}}$ используются для вычисления компонент $z_{3}^{{ijk}}$. При этом набор индексов $k$, которые используются для вычисления компонент $z_{1}^{{ijk}}$, зависит от конкретного индекса $j$ компоненты $z_{3}^{{ijk}}$. Из этих соображений мы присваиваем компонентам $z_{1}^{{ijk}}$ и индекс $j$, хотя явно они от него не зависят. Аналогично, и компоненты $z_{2}^{{ijk}}$ не зависят от индекса $j$, но мы им присваиваем этот индекс.

Далее для всех фиксированных пар индексов $(i,j)$, таких что $i = \overline {1,I} $_, $j = \overline {1,I} $ и $j \ne i$, вычисляем компоненты $z_{9}^{{ij}},$ $z_{{10}}^{{ij}},$ $z_{{11}}^{{ij}},$ $z_{{12}}^{{ij}},$ $z_{{14}}^{{ij}},$ $z_{{15}}^{{ij}},$ $z_{{16}}^{{ij}},$ $z_{{17}}^{{ij}}$.Компоненты $z_{{13}}^{{ij}}$ также не зависят от третьего индекса ($k$), но они используются для вычисления компонент $z_{3}^{{ijk}}$, поэтому вычисляются раньше.

Наконец, энергия системы атомов вычисляется по формуле

где фазовые переменные $z_{1}^{{}},z_{2}^{{}},...,z_{{17}}^{{}}$ определяются указанным выше многошаговым процессом.

Для получения аналитических выражений производных функции $E$ по координатам атомов воспользуемся структурой функции $E$, которая представлена формулой (1.10). Для всех фиксированных пар индексов $(l,m)$, таких что $l = 1,\;2,\;3$ и $m = \overline {1,I} $, введем следующие обозначения: $\tilde {z}_{1}^{{}},\tilde {z}_{2}^{{}},...,\tilde {z}_{{17}}^{{}}$, где

Указанные выше величины вычисляются по следующим формулам:

– для всех $l = 1,\;2,\;3$ и $i,j,k,m = \overline {1,I} $

– для всех $l = 1,\;2,\;3$ и $i,j,m = \overline {1,I} $

Порядок вычисления этих производных такой же, как и для фазовых переменных $z_{1}^{{}},z_{2}^{{}},...,z_{{17}}^{{}}$, а формулы для вычисления частных производных функции энергии по координатам атомов имеют вид

$\frac{{\partial E(u,x)}}{{\partial {{x}_{{lm}}}}} = \sum\nolimits_{i = 1}^I {\sum\nolimits_{\begin{subarray}{l} j = 1 \\ j \ne i \end{subarray}} ^I {\tilde {z}_{{17}}^{{ij}}} } $,  $l = 1,\;2,\;3$  и  $m = \overline {1,I} $.

Проведена большая серия вычислительных экспериментов, в которых менялись размерность материала (рассматривались трехмерная и двумерная модели) и количество атомов в фрагменте материала. Если рассматривается трехмерная модель $(l = 1,2,3)$, то для каждого $m$-го атома вычисления по многошаговому процессу $\tilde {z}_{1}^{{}},\tilde {z}_{2}^{{}},...,\tilde {z}_{{17}}^{{}}$ необходимо проводить 3 раза, а если рассматривается двумерная модель, то 2 раза. Все численные эксперименты показали, что в первом случае (трехмерная модель) градиент функции энергии по координатам атома, определенный с помощью формул БАД-методологии (см. [12]), вычисляется в 3 раза быстрее, чем с помощью многошагового процесса $\tilde {z}_{1}^{{}},\tilde {z}_{2}^{{}},...,\tilde {z}_{{17}}^{{}}$. Когда используется двумерная модель материала, эффективность формул БАД-методологии возрастает в 2 раза.

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

В [13] рассмотрена задача идентификации параметров потенциала Терсоффа применительно к однокомпонентному кристаллу кремния. Для определения начальных приближений в указанной работе допустимое множество $X = [\underline u ,\bar {u}] = \left\{ {u \in {{R}^{{10}}}:{{{\underline u }}_{i}} \leqslant {{x}_{i}} \leqslant {{{\bar {u}}}_{i}}} \right\}$ выбиралось таким образом, чтобы заведомо содержались все допустимые значения параметров, а именно:

В табл. 1–3 представлены результаты сравнения вычисления градиента потенциала Терсоффа на основе формул БАД-методологии и вычисления этого градиента с помощью метода конечных разностей для трех компонент градиента (по параметрам ${{u}_{1}} = {{D}_{e}},$ ${{u}_{3}} = \beta ,$ ${{u}_{6}} = \gamma $). Все входные параметры задачи выбирались из указанного выше диапазона допустимых значений. Величина $\Delta $ в табл. 1–3 указывает на значения приращения параметра ${{u}_{i}}$, которые были использованы при вычислении градиента с помощью метода конечных разностей.

В табл. 1 сравниваются значения первой компоненты градиента энергии по параметру ${{u}_{1}} = {{D}_{e}}$ при ${{D}_{e}} = 1$. Значение этой компоненты градиента, вычисленное с помощью БАД-методологии (точное значение), равно $\frac{{\partial E(u)}}{{\partial {{u}_{1}}}} = - {\text{7}}{\text{.079590296064}}$.

В табл. 2 сравниваются значения третьей компоненты градиента энергии по параметру ${{u}_{3}} = \beta $ при $\beta = 1$. Значение этой компоненты градиента, вычисленное с помощью БАД-методологии (точное значение), равно $\frac{{\partial E(u)}}{{\partial {{u}_{3}}}} = {\text{8}}{\text{.117021018694}}$.

В табл. 3 сравниваются значения шестой компоненты градиента энергии по параметру ${{u}_{6}} = \gamma $ при γ = 1.0 × 10^–7. Значение этой компоненты градиента, вычисленное с помощью БАД-методологии (точное значение), равно $\frac{{\partial E(u)}}{{\partial {{u}_{6}}}} = {\text{1}}{\text{.288320939087}} \times {\text{1}}{{{\text{0}}}^{{ + {\text{6}}}}}$.

Как следует из табл. 1, при вычислении первой компоненты градиента наиболее подходящим является выбор $\Delta $ из интервала $[{{10}^{{ - 8}}};\;1.0]$, что составляет от 10^–7 до $100\% $ от начального приближения ${{u}_{1}} = 1.0$.

Из табл. 2 следует, что при вычислении третьей компоненты градиента наиболее подходящим является выбор $\Delta \approx {{10}^{{ - 7}}}$. Это составляет примерно ${{10}^{{ - 5}}}\% $ от начального приближения ${{u}_{3}} = 1.0$.

Из табл. 3 следует, что при вычислении шестой компоненты градиента наиболее подходящим является выбор $\Delta \approx {{10}^{{ - 11}}}$. Это составляет примерно ${{10}^{{ - 2}}}\% $ от начального приближения ${{u}_{6}} = {{10}^{{ - 7}}}$.

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

2. СРАВНЕНИЕ МЕТОДОВ ОПРЕДЕЛЕНИЯ ПРОИЗВОДНЫХ ДЛЯ ПЛОСКОЙ МОДЕЛИ МНОГОСЛОЙНОГО КУСОЧНО-ОДНОРОДНОГО МАТЕРИАЛА

В качестве вычислительного эксперимента для сравнения методов определения производных с помощью указанных выше способов и с помощью стандартных пакетов программ рассмотрена предложенная в [14] плоская модель многослойного кусочно-однородного материала. Материал представляется в виде периодической кусочно-однородной многослойной структуры, в рамках которой могут отличаться типы атомов в различных слоях. Модель накладывает ряд ограничений на структуру слоев:

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

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

– в рассматриваемой системе выделяется группа из $K$ параллельных слоев, которая периодически повторяется в направлении оси $y$ (см. фиг. 1);

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

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

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

${{h}_{i}},\;\,\,i = 1,...,K$, – расстояние между слоем с номером $i$ и предыдущим слоем,

${{d}_{i}},\,\;\,i = 1,...,K$, – смещение первого атома в слое $i$ с положительной абсциссой относительно нулевой отметки,

${{s}_{i}},\,\;\,i = 1,...,K$, – расстояние между атомами в слое $i$.

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

Оптимизационная задача состоит в минимизации энергии $E(u,w)$ группы атомов (потенциал взаимодействия атомов – потенциал Терсоффа), расположенных на $K$ соседних слоях. Параметрами оптимизационной задачи являются переменные $w = \left( {{{h}_{1}},{{d}_{1}},{{s}_{1}},...,{{h}_{K}},{{d}_{K}},{{s}_{K}}} \right)$, составляющие вектор размерности $r = 3K$.

Пусть (для наглядности) ${{x}_{i}} = {{x}_{{1i}}}$ – первая координата $i$-го атома двумерной структуры, ${{y}_{i}} = {{x}_{{2i}}}$ – его вторая координата. Если $i$-й атом является атомом с порядковым номером $j$ на $k$‑м слое, то

Соответствующие производные вычисляются по формулам (см. разд. 1) где $\frac{{\partial E}}{{\partial {{x}_{i}}}}$ и $\frac{{\partial E}}{{\partial {{y}_{i}}}}$ – производные функции энергии по первой и второй координатам атома соответственно.

Расчет производных энергии $E(u,w)$ по параметрам $w = \left( {{{h}_{1}},{{d}_{1}},{{s}_{1}},...,{{h}_{K}},{{d}_{K}},{{s}_{K}}} \right)$ осуществлялся как с использованием аналитически выведенных формул, так и с помощью имеющихся пакетов программ.

Работа выполнялась с использованием инфраструктуры Центра коллективного пользования “Высокопроизводительные вычисления и большие данные” (ЦКП “Информатика”) ФИЦ ИУ РАН (г. Москва). Испытательный стенд состоял из следующих средств вычислительной техники:

• Intel Core i7 4770 3.4 Гц, 8 Гб RAM – 4-х ядерный процессор архитектуры Haswell c поддержкой технологий Hyper-Threading (два вычислительных протока на ядро, всего восемь логических ядер), AVX2 –Advanced Vector Extensions (в контексте решаемой задачи – обработка данных в формате с плавающей запятой в группах длиной 256 бит, т.е. для чисел с двойной точностью одновременно обрабатывается вектор из четырех чисел), и FMA–Fused Multiply-Add – выполнение совмещенной операции умножения-сложения вида a = a + b*c. Использовался компилятор GNUC++ 6.4.0.

• 2x Intel Xeon E5-2683 V4 2.1 Гц, 512 Гб RAM – двухпроцессорный сервер, каждый процессор содержит 16 ядер с поддержкой Hyper-Threading, AVX2 и FMA (всего 64 логических ядра). Использовался компилятор GNUC++ 7.2.0.

В работе использовались хорошо зарекомендовавшие себя ранее (см. [15]) пакеты CoDiPack 1.6.0 (см. [10]) и Adept 1.1 (см. [11]).

Сравнивались следующие методы вычисления градиента сложной функции $E(u,w)$:

1) использование выведенных аналитически выражений для вычисления частных производных этой функции (метод Hand1),

2) использование выведенных аналитически формул БАД-методологии для вычисления частных производных этой функции (метод Hand2),

3) использование прямого метода автоматического дифференцирования, реализованного в пакете CoDiPack с типом данных RealForwardVecGen (метод CoDiPack_kr),

4) использование прямого метода автоматического дифференцирования, реализованного в пакете CoDiPack с типом данных RealForward (метод CoDiPack_rk),

5) использование прямого метода автоматического дифференцирования, реализованного в пакете CoDiPack с типом данных RealForward, распараллеленный с помощью директивы OpenMP ‑ #pragma omp for (метод CoDiPack_rk),

6) использование обратного метода автоматического дифференцирования, реализованного в пакете Adept ([11]) (метод Adept).

Для методов Hand1, CoDiPack_kr, CoDiPack_rk и CoDiPack_rk_omp использовались опции оптимизации компилятора O2 и Ofast. Для методов Hand2 и Adept только опция O2.

Расчеты проводились для фиксированного набора специфических параметров потенциала Терсоффа, но для различного количества слоев и атомов:

– 4 слоя и 94 атома, размерность $r$ вектора $w$ равна 12 (распределение атомов по слоям 28, 19, 28, 19),

– 8 слоев и 188 атомов, размерность $r$ вектора $w$ равна 24 (распределение атомов по слоям 28, 19, 19, 28, 28, 19, 19, 28),

– 12 слоев и 204 атома, размерность $r$ вектора $w$ равна 36 (распределение атомов по слоям 15, 19, 19, 15, 15, 19, 19, 15, 15, 19, 19, 15).

При расчете градиента пакетом CoDiPack применялись типы данных RealForward и RealForwardVecGen. Тип данных RealForward состоит из двух чисел – значения переменной и ее производной по выбранной переменной. Для расчета градиента функции необходимо запустить вычисление функции и производной в цикле $r$ раз. Тип данных RealForwardVecGen состоит из $r + 1$ чисел – значения переменной и вектора ее частных производных. При расчете градиента цикл по $r$ становится внутренним, но функция рассчитывается один раз. Использование типа данных RealForwardVecGen позволяет задействовать технологию AVX2 и ускоряет вычисления от 1.5 до 4 раз без каких-либо трудозатрат; достаточно указать компилятору опцию оптимизации –Ofast.

В табл. 4 приводится сравнение прямых методов дифференцирования Hand1, CoDiPack_kr, CoDiPack_rk в однопоточном исполнении и CoDiPack_rk_omp, распараллеленный с помощью директивы OpenMP – #pragma omp for. Как видно из таблицы, время вычисления функции имеет примерно квадратичную зависимость от количества атомов.

Время вычисления производной методом Hand1 изменяется на 10–17% в зависимости от установленных опций оптимизации, при этом для IntelCore i7 опция Ofast ускоряет расчет, а для IntelXeon наоборот замедляет.

Метод CoDiPack_kr при опции O2 оказывается в 3–4 раза быстрее метода Hand1 с той же опцией. Связано это с тем, что при реализации методов Hand1 и Hand2 задача оптимизации кода не ставилась, а пакетные методы оптимизировались. Сравнивая влияние опций O2 и Ofast, мы видим существенное ускорение решения задачи при $r = 12$ и несколько уменьшающееся ускорение для больших размерностей. Для малой размерности задачи преимущества технологии AVX2 раскрываются полностью.

Метод CoDiPack_rk, как можно увидеть, не так сильно чувствителен к опциям оптимизации компилятора. Вычислительно он более сложен, так как значение функции рассчитывается $r$ раз. Но его применение целесообразно при распараллеливании процесса расчета градиента на несколько ядер. Эксперименты с методом CoDiPack_rk_omp это показывают.

В табл. 5 приведено время вычисления градиента и функции обратными методами Hand2 и Adept. Как видно из таблицы, время вычисления функции и градиента мало зависит от размерности задачи. Исключение составляет Hand2 для задачи с $r = 12$ (4*94). Связано это с небольшим количеством памяти, при котором вычисления выполняются в основном с использованием кэшей процессора. Как и в случае с прямыми методами оптимизация кода Hand2 не проводилась, с чем и связано большее время вычислений. Также необходимо обратить внимание на тот факт, что вычисления прямым методом CoDiPack_rk_omp на компьютере IntelXeon оказываются более эффективными, чем вычисления обратным методом Adept, за счет параллельного исполнения.

3. ЗАКЛЮЧЕНИЕ

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

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

Что касается “ручной” версии БАД-методологии, то в этом случае необходимо вручную построить многошаговый процесс вычисления функции, с помощью канонических формул БАД‑методологии построить многошаговый процесс решения сопряженной задачи и вывести формулы для вычисления градиента функции. Безусловно, все это требует серьезных интеллектуальных затрат исследователя. Однако при этом у исследователя появляется возможность учесть специфику конкретной задачи и сделать алгоритм вычисления градиента функции более эффективным. К преимуществам “ручной” версии БАД-методологии следует отнести и то, что эта версия гораздо слабее, чем “пакетная” версия, зависит от параметров используемого компьютера.