Сравнение гладких и негладких методов решения задачи о нахождении пространственной структуры молекулы

Литвинцев Игорь Сергеевич

Сравнение гладких и негладких методов решения задачи о нахождении пространственной структуры молекулы

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

Общественные науки в целом

Дипломы

Вуз: Санкт-Петербургский государственный университет (СПбГУ)

ID: 587d364b5f1be77c40d58bc7

UUID: 2d8b8038-7322-4574-b108-23d70e2ec170

Язык: Русский

Опубликовано: почти 8 лет назад

Просмотры: 8

Литвинцев Игорь Сергеевич

Источник: Санкт-Петербургский государственный университет

Комментировать 0

Рецензировать 0

Скачать - 425628 bytes

Поделиться работой

САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ КАФЕДРА МАТЕМАТИЧЕСКОЙ ТЕОРИИ МОДЕЛИРОВАНИЯ СИСТЕМ УПРАВЛЕНИЯ Литвинцев Игорь Сергеевич Выпускная квалификационная работа бакалавра СРАВНЕНИЕ ГЛАДКИХ И НЕГЛАДКИХ МЕТОДОВ РЕШЕНИЯ ЗАДАЧИ О НАХОЖДЕНИИ ПРОСТРАНСТВЕННОЙ СТРУКТУРЫ МОЛЕКУЛЫ Направление 010400 Прикладная математика, фундаментальная информатика и основы программирования . Научный руководитель, кандидат физ.-мат. наук, доцент Тамасян Г. Ш. Санкт-Петербург 2016

Содержание Введение 3 1 Постановка задачи 4 2 Вспомогательные сведения 5 2.1 Определения и обозначения . . . . . . . . . . . . . . . . . . . 5 2.2 Элементы гиподифференциального исчисления . . . . . . . . 6 2.3 МДМ-метод . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.4 Алгоритмы минимизации . . . . . . . . . . . . . . . . . . . . 8 3 4 Решение задачи 11 3.1 Минимизация гладкого функционала . . . . . . . . . . . . . 11 3.2 Минимизация негладких функционалов . . . . . . . . . . . . 13 Численные эксперименты 18 Заключение 22 Список литературы 23 Приложение 24 2

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

1 Постановка задачи Задача о нахождении пространственной структуры молекулы может быть сформулирована, как задача поиска декартовых координат атомов 𝑋1 , . . . , 𝑋𝑚 молекулы таких, что ‖𝑋𝑖 − 𝑋𝑗 ‖ = 𝑑𝑖𝑗 ([𝑖, 𝑗] ∈ 𝐷), (1) где 𝐷 набор пар атомов [𝑖, 𝑗], между которыми известно евклидово расстояние 𝑑𝑖𝑗 . Множество 𝐷 можно представить, как симметрическую матрицу, в которой элемент 𝑑𝑖𝑗 — расстояние между 𝑖-ым и 𝑗 -ым атомами. Если расстояние между некоторой парой атомов неизвестно, на соответствующей позиции в матрице 𝐷 поставим ноль. В случае если точные расстояния между всеми парами атомов заданы, решение задачи, согласно статье [1], находится за линейное время. Однако на практике зачастую не только не все межатомные расстояния известны, но и не всегда существует возможность измерить расстояние между некоторыми парами атомов. Есть множество различных подходов к решению задачи с разреженной матрицей расстояний. Один из них основан на том, что корректная конфигурация молекулы должна совпадать с минимумом некоторой функции невязок. Введем три различные функции невязок ∑︁ (‖𝑋𝑖 − 𝑋𝑗 ‖2 − 𝑑2𝑖𝑗 )2 , (2) |‖𝑋𝑖 − 𝑋𝑗 ‖2 − 𝑑2𝑖𝑗 |, (3) 𝐹3 (𝑋1 , . . . , 𝑋𝑚 ) = max |‖𝑋𝑖 − 𝑋𝑗 ‖2 − 𝑑2𝑖𝑗 |, (4) 𝐹1 (𝑋1 , . . . , 𝑋𝑚 ) = [𝑖,𝑗]∈𝐷 ∑︁ 𝐹2 (𝑋1 , . . . , 𝑋𝑚 ) = [𝑖,𝑗]∈𝐷 [𝑖,𝑗]∈𝐷 где 𝑚 – число атомов в молекуле. Заметим, что первая функция является гладкой, а две другие нет. Реализация и сравнение алгоритмов минимизации приведенных функционалов составляют цель данной работы. 4

2 Вспомогательные сведения Все нижеизложенные сведения данной главы взяты из [4]–[8]. 2.1 Определения и обозначения Определение 2.1. Выпуклой оболочкой точек 𝑋1 , . . . , 𝑋𝑚 называется мно- жество {︃ co {𝑋1 , . . . , 𝑋𝑚 } = 𝑋= 𝑚 ∑︁ 𝑘=1 Определение 2.2. ⃒ }︃ 𝑚 ⃒ ∑︁ ⃒ 𝛼𝑘 = 1 . 𝛼𝑘 𝑋𝑘 ⃒ 𝛼𝑘 > 0, ⃒ 𝑘=1 Функция 𝑓 , заданная и конечная на открытом множе- стве 𝑋 ⊂ R𝑛 , гиподифференцируема в точке 𝑥 ∈ 𝑋 , если существует такой выпуклый компакт 𝑑𝑓 (𝑥) ⊂ R𝑛+1 , что 𝑓 (𝑥 + Δ) = 𝑓 (𝑥) + max [𝑎 + ⟨𝑣, Δ⟩] + 𝑜𝑥 (Δ), [𝑎,𝑣]∈𝑑𝑓 (𝑥) где lim 𝛼↓0 𝑜𝑥 (𝛼, Δ) =0 𝛼 ∀Δ ∈ R𝑛 . Здесь 𝑎 ∈ R1 ; 𝑣 ∈ R𝑛 . Множество 𝑑𝑓 (𝑥) называется гиподифференциалом функции 𝑓 в точке 𝑥. Определение 2.3. Пусть 𝐴, 𝐵 ⊂ R𝑛 . Число 𝜌(𝐴, 𝐵) = sup inf ‖𝑣 − 𝑤‖ + sup inf ‖𝑣 − 𝑤‖ 𝑤∈𝐴 𝑣∈𝐵 𝑣∈𝐵 𝑤∈𝐴 называется расстоянием между множествами в метрике Хаусдорфа. Определение 2.4. Отображение 𝐴(𝑥) называется непрерывным по Хау- сдорфу в точке 𝑥0 , если lim 𝜌(𝐴(𝑥), 𝐴(𝑥0 )) = 0 𝑥→𝑥0 Определение 2.5. Функция 𝑓 называется непрерывно гиподифференци- руемой в точке 𝑥, если она гиподифференцируема в некоторой окрестности точки 𝑥 и если существует непрерывное (по Хаусдорфу) в этой точке гиподифференциальное отображение 𝑑𝑓 . 5

2.2 Элементы гиподифференциального исчисления Везде ниже все рассматриваемые функции полагаются заданными и конечными на открытом множестве 𝑋 ⊂ R𝑛 . А суммирование выпуклых компактов производится по Минковскому. Лемма 2.1. Гиподифференциал гладкой функции 𝑓 (𝑥) в точке 𝑥 ∈ 𝑋 равен {︃(︃ 𝑑𝑓 (𝑥) = Лемма 2.2. 0 )︃}︃ ∇𝑓 (𝑥) . (5) Если 𝑓1 , . . . , 𝑓𝑚 — гиподифференцируемые в точке 𝑥 ∈ 𝑋 функции, то и функция 𝑓 = ема в этой точке, при этом 𝑚 ∑︀ 𝑐𝑖 𝑓𝑖 , где 𝑐𝑖 ∈ R1+ , тоже гиподифференциру- 𝑖=1 𝑑𝑓 (𝑥) = 𝑚 ∑︁ 𝑐𝑖 𝑑𝑓𝑖 (𝑥). (6) 𝑖=1 Лемма 2.3. Если 𝑓𝑖 , при 𝑖 ∈ 1 : 𝑚 — гиподифференцируемые в точке 𝑥 функции, то и функция 𝑓 (𝑥) = max 𝑓𝑖 (𝑥) тоже гиподифференцируема в 𝑖∈1:𝑚 этой точке, при этом {︃ (︃ 𝑑𝑓 (𝑥) = co 𝑑𝑓𝑖 (𝑥) + 2.3 )︃ ⃒ }︃ ⃒ 𝑓𝑖 (𝑥) − 𝑓 (𝑥) ⃒ ⃒ 𝑖 ∈ 1:𝑚 . ⃒ 0𝑛 (7) МДМ-метод Пусть в пространстве R𝑛 заданы 𝑚 точек, 𝐻 = {𝑎𝑖 }𝑚 𝑖=1 . Обозначим через 𝐺 выпуклую оболочку множества 𝐻 . Задача: найти точку из 𝐺, ближайшую к началу координат или, что эквивалентно, решить следующую задачу ‖𝑣‖2 → min . 𝑣∈𝐺 (8) В силу единственности точки минимума выпуклой функции на выпуклом множестве, задача (8) имеет единственное решение. Обозначим его 𝑣* . 6

Обозначим через 𝐴 матрицу со столбцами 𝑎1 , . . . , 𝑎𝑚 . Тогда любой вектор 𝑣 ∈ 𝐺 допускает представление 𝑣 = 𝐴𝑝, 𝑝 ≥ 0𝑚 , 𝑚 ∑︁ 𝑝(𝑖) = 1. (9) 𝑖=1 Обозначим 𝑀+ (𝑝) носитель вектора 𝑝 𝑀+ (𝑝) = {𝑖 ∈ 1 : 𝑚 | 𝑝(𝑖) > 0}. Введем величину Δ(𝑝) = max ⟨𝑎𝑖 , 𝑣⟩ − min ⟨𝑎𝑖 , 𝑣⟩ , 𝑖∈1:𝑚 𝑖∈𝑀+ (𝑝) где 𝑣 = 𝐴𝑝. Вектор 𝑝 удовлетворяет условиям в (9). Множество таких векторов обозначим 𝑃 . Алгоритм 1) Выберем любое 𝑣0 ∈ R𝑛 . 2) Пусть уже построено 𝑣𝑘 = 𝐴𝑝𝑘 , 𝑝𝑘 ∈ 𝑃 . Найдем индексы 𝑖′𝑘 ∈ 𝑀+ (𝑝𝑘 ) и 𝑖′′𝑘 ∈ 1 : 𝑚, такие, что ⟨︀ ⟩︀ max ⟨𝑎𝑖 , 𝑣𝑘 ⟩ = 𝑎𝑖′𝑘 , 𝑣𝑘 , 𝑖∈𝑀+ (𝑝𝑘 ) ⟨︀ ⟩︀ min ⟨𝑎𝑖 , 𝑣𝑘 ⟩ = 𝑎𝑖′′𝑘 , 𝑣𝑘 . 𝑖∈1:𝑚 Обозначим 𝑎𝑖′𝑘 = 𝑎′𝑘 , 𝑎𝑖′′𝑘 = 𝑎′′𝑘 . Найдем Δ𝑘 := Δ(𝑝𝑘 ) = ⟨𝑎′𝑘 − 𝑎′′𝑘 , 𝑣𝑘 ⟩ . Если Δ𝑘 = 0, то 𝑣𝑘 — решение задачи (8). Процесс закончен. 3) Пусть Δ𝑘 > 0. Следующее приближение ищем, минимизируя выраже- ние 𝑣𝑘 (𝑡) = 𝑣𝑘 − 𝑡𝑝′𝑘 (𝑎′𝑘 − 𝑎′′𝑘 ), 𝑡 ∈ [0, 1], (𝑖′ ) где 𝑝′𝑘 = 𝑝𝑘 𝑘 . Выберем 𝑡𝑘 из условия ‖𝑣𝑘 (𝑡𝑘 )‖2 = min ‖𝑣𝑘 (𝑡)‖2 . 𝑡∈[0,1] 7

Укажем явную формулу для 𝑡𝑘 ⎧ ⎨̂︀ 𝑡𝑘 , при ̂︀ 𝑡𝑘 < 1; 𝑡𝑘 = ⎩1, при ̂︀ 𝑡𝑘 ≥ 1; где Δ𝑘 𝑝′𝑘 Δ𝑘 = . ′ ′ ‖̂︀ 𝑣𝑘 − 𝑣𝑘 ‖2 𝑝𝑘 ‖𝑎𝑘 − 𝑎′′𝑘 ‖2 Вектор 𝑝𝑘+1 находим по формуле ⎧ ⎪ (𝑖) ⎪ 𝑝 при 𝑖 ̸= 𝑖′𝑘 , 𝑖 ̸= 𝑖′′𝑘 , ⎪ 𝑘 , ⎪ ⎨ ′ (𝑖) 𝑘) 𝑝𝑘+1 = (1 − 𝑡𝑘 )𝑝(𝑖 , при 𝑖 = 𝑖′𝑘 , 𝑘 ⎪ ⎪ ⎪ ′ ′′ ⎪ ⎩𝑝(𝑖𝑘 ) + 𝑡𝑘 𝑝(𝑖𝑘 ) , при 𝑖 = 𝑖′′ . 𝑘 𝑘 𝑘 ̂︀ 𝑡𝑘 = 4) 5) Положим 𝑣𝑘+1 = 𝑣𝑘 (𝑡𝑘 ). Описание МДМ-метода завершено. Последовательность {𝑣𝑘 } сходит- ся к 𝑣* — решению задачи (8). 2.4 Алгоритмы минимизации Метод градиентного спуска Пусть функция 𝑓 определена и дифференцируема в пространстве R𝑛 . Будем осуществлять минимизацию в направлении наискорейшего спуска, т.е. в направлении антиградиента −∇𝑓 . Алгоритм 1) Выберем любое 𝑥0 ∈ R𝑛 . 2) Пусть уже построено 𝑥𝑘 ∈ R𝑛 . Тогда 𝑥𝑘+1 = 𝑥𝑘 − 𝛼𝑘 ∇𝑓 (𝑥𝑘 ), (10) где 𝛼𝑘 — находим, решая задачу одномерной минимизации 𝛼𝑘 = argmin 𝑓 (𝑥𝑘 − 𝛼∇𝑓 (𝑥𝑘 )). (11) 𝛼∈[0,∞] 3) Если для заданного 𝜀 выполнено одно из условий ‖𝑥𝑘+1 − 𝑥𝑘 ‖ < 𝜀; (12) ‖𝑓 (𝑥𝑘+1 ) − 𝑓 (𝑥𝑘 )‖ < 𝜀; (13) ‖∇𝑓 (𝑥𝑘+1 )‖ < 𝜀. (14) 8

то 𝑥* = 𝑥𝑘+1 — точка минимума и алгоритм завершен. Иначе переходим к следующей итерации. Метод гиподифференциального спуска Пусть функция 𝑓 задана, липшицева и непрерывно гиподифференцируема на R𝑛 . Лемма 2.4. Необходимое условие минимума гиподифференцируемой функ- ции 𝑓 в точке 𝑥 имеет вид 0𝑛+1 ∈ 𝑑𝑓 (𝑥). (15) Алгоритм 1) Выберем любое 𝑥0 ∈ R𝑛 . 2) Пусть уже построено 𝑥𝑘 ∈ R𝑛 . Если в точке 𝑥𝑘 выполнено (15), то 𝑥𝑘 — точка минимума, и процесс останавливается. 3) Иначе найдем min ‖𝑧‖ = 𝑧¯𝑘 . 𝑧∈𝑑𝑓 (𝑥𝑘 ) (16) Из того, что 𝑥𝑘 не является точкой минимума следует, что 𝑧¯𝑘 = [𝜂𝑘 , 𝑧𝑘 ] ̸= 0𝑛+1 , где 𝜂𝑘 ∈ R1 , 𝑧𝑘 ∈ R𝑛 . Направление наискорейшего спуска находится по формуле 𝑔𝑘 = − 4) 𝑧𝑘 . ‖𝑧𝑘 ‖ Решая задачу одномерной минимизации, найдем шаг 𝛼𝑘 𝛼𝑘 = argmin 𝑓 (𝑥𝑘 + 𝛼𝑔𝑘 ). 𝛼>0 5) (17) Следующую точку положим 𝑥𝑘+1 = 𝑥𝑘 + 𝛼𝑘 𝑔𝑘 . Далее продолжаем аналогично. В результате получаем последовательность {𝑥𝑘 } такую, что 𝑓 (𝑥𝑘+1 ) < 𝑓 (𝑥𝑘 ). 9

На практике условие (15) будем проверять следующим образом. Если для заданного 𝜀 выполнено одно из условий ‖𝑥𝑘+1 − 𝑥𝑘 ‖ < 𝜀; (18) ‖𝑓 (𝑥𝑘+1 ) − 𝑓 (𝑥𝑘 )‖ < 𝜀; (19) ‖¯ 𝑧𝑘 ‖ < 𝜀; (20) то 𝑥* = 𝑥𝑘+1 — точка минимума и алгоритм завершен. Иначе переходим к следующей итерации. 10

3 Решение задачи Обозначим координаты точки 𝑋𝑖 ∈ R3 , как (𝑥𝑖 , 𝑦𝑖 , 𝑧𝑖 ), где 𝑖 ∈ 1 : 𝑚 (случай 𝑋𝑖 ∈ R2 аналогичен). Введем вектор переменных 𝑥 = (𝑋1 , . . . , 𝑋𝑚 ) = (𝑥1 , 𝑦1 , 𝑧1 , . . . , 𝑥𝑚 , 𝑦𝑚 , 𝑧𝑚 ). Без ограничения общности первый атом можно связать с началом координат. Таким образом, в дальнейшем при реализации алгоритмов оптимизации будем исключать координаты первого атома 𝑋𝑚 = (0, 0, 0) из искомых переменных. Следовательно, количество переменных равняется 𝑁 = 𝑛(𝑚 − 1), где 𝑛 – размерность пространства (везде далее будем полагать 𝑛 = 3). 3.1 Минимизация гладкого функционала Для минимизации функционала (2) применим метод наискорейшего спуска. Выпишем формулу, по которой будем вычислять градиент в заданной точке. Для начала перепишем функционал в векторном виде 𝐹1 (𝑋1 , . . . , 𝑋𝑚 ) = ∑︁ (︀ )︀2 (𝑋𝑖 − 𝑋𝑗 )𝑇 (𝑋𝑖 − 𝑋𝑗 ) − 𝑑2𝑖𝑗 , 1≤𝑖<𝑗≤𝑚 𝑑𝑖𝑗 ̸=0 Перегруппируем слагаемые следующим образом 𝐹1 (𝑋1 , . . . , 𝑋𝑚 ) = 𝑚 ∑︁ (︀ )︀2 (𝑋1 − 𝑋𝑗 )𝑇 (𝑋1 − 𝑋𝑗 ) − 𝑑21𝑗 + 𝑗=2 𝑑1𝑗 ̸=0 + 𝑚 ∑︁ (︀ )︀2 (𝑋2 − 𝑋𝑗 )𝑇 (𝑋2 − 𝑋𝑗 ) − 𝑑22𝑗 + · · · + 𝑗=3 𝑑2𝑗 ̸=0 + 𝑚 ∑︁ (︀ )︀2 (𝑋𝑚−2 − 𝑋𝑗 )𝑇 (𝑋𝑚−2 − 𝑋𝑗 ) − 𝑑2𝑚−2,𝑗 + 𝑗=𝑚−1 𝑑𝑚−2,𝑗 ̸=0 + 𝑚 ∑︁ (︀ )︀2 (𝑋𝑚−1 − 𝑋𝑗 )𝑇 (𝑋𝑚−1 − 𝑋𝑗 ) − 𝑑2𝑚−1,𝑗 . 𝑗=𝑚 𝑑𝑚−1,𝑗 ̸=0 11

Обозначим 𝜑𝑖 (𝑋𝑖 , . . . , 𝑋𝑚 ) = 𝑚 ∑︁ (︀ )︀2 (𝑋𝑖 − 𝑋𝑗 )𝑇 (𝑋𝑖 − 𝑋𝑗 ) − 𝑑2𝑖𝑗 , ∀𝑖 ∈ 1 : 𝑚−1. 𝑗=𝑖+1 𝑑𝑖𝑗 ̸=0 Отметим, что каждое 𝜑𝑖 при 𝑖 ∈ 1 : 𝑚 − 1 зависит только от точек 𝑋𝑖 , . . . , 𝑋 𝑚 . Обозначим за 𝜕 𝜕𝑋𝑘 вектор частных производных по переменным 𝑥𝑘 , 𝑦𝑘 , 𝑧𝑘 𝜕 = 𝜕𝑋𝑘 Тогда для 𝜕𝜑𝑖 𝜕𝑋𝑘 (︂ 𝜕 𝜕 𝜕 , , 𝜕𝑥𝑘 𝜕𝑦𝑘 𝜕𝑧𝑘 )︂ . справедливы формулы 𝜕𝜑𝑖 = 0𝑛 , 𝜕𝑋𝑘 ∀𝑖 ∈ 2 : 𝑚−1, 𝑘 ∈ 1 : 𝑖−1, 𝑚 ∑︁ (︀ )︀ 𝜕𝜑𝑖 (𝑋𝑖 − 𝑋𝑗 )𝑇 (𝑋𝑖 − 𝑋𝑗 ) − 𝑑2𝑖𝑗 (𝑋𝑖 − 𝑋𝑗 ), =4 𝜕𝑋𝑖 𝑗=𝑖+1 ∀𝑖 ∈ 1 : 𝑚−1, 𝑑𝑖𝑗 ̸=0 (︀ )︀ 𝜕𝜑𝑖 = −4 (𝑋𝑖 − 𝑋𝑘 )𝑇 (𝑋𝑖 − 𝑋𝑘 ) − 𝑑2𝑖𝑘 (𝑋𝑖 − 𝑋𝑘 ), 𝜕𝑋𝑘 ∀𝑖 ∈ 1 : 𝑚−1, 𝑘 ∈ 𝑖+1 : 𝑚, 𝑑𝑖𝑘 ̸= 0. Наконец, компоненты градиента функции 𝐹1 (𝑋1 , . . . , 𝑋𝑚 ) находятся по следующей формуле 𝑚−1 ∑︁ 𝜕𝜑𝑖 𝜕𝐹1 = 𝜕𝑋𝑘 𝜕𝑋𝑘 𝑖=1 ∀𝑘 ∈ 1 : 𝑚. Учитывая, что точку 𝑋1 мы зафиксировали в начале координат, в получившемся градиенте отбросим первые 𝑛 компонент, соответствующие частным производным по компонентам точки 𝑋1 . Таким образом, ∇𝐹1 будет иметь размерность 𝑁 = 𝑛(𝑚−1). Задачу одномерной минимизации (11) будем решать, используя встроенную в MATLAB функцию fmincon, реализующую алгоритм внутренней точки. Отрезок минимизации формируем следующим образом: на первой итерации он равен фиксированному отрезку [0, 10𝑑𝑚𝑎𝑥 ], где 𝑑𝑚𝑎𝑥 — максимальный элемент матрицы расстояний 𝐷. На 𝑘 -й итерации (𝑘 > 1) отрезок 12

* * становится равен [0, 5𝛼𝑘−1 ], где 𝛼𝑘−1 — минимум задачи (11) на (𝑘 −1)-й итерации. В качестве правила останова используем каждый из критериев (12)– (14) по очереди и сведем в единую таблицу количество итераций прошедших до выхода. 3.2 Минимизация негладких функционалов Для минимизации негладких функционалов (3) и (4) будем использовать метод гиподифференциального спуска. Минимизация функционала 𝐹2 (𝑥) Для начала опишем процесс вычисления гиподифференциала функции 𝐹2 в произвольно заданной точке 𝑥. Перепишем модули под знаком суммы в (3), как максимум двух функций 𝐹2 (𝑥) = 𝑀 ∑︁ 𝑓𝑘𝑖𝑗 (𝑥) = 𝑘𝑖𝑗 =1 𝑀 ∑︁ max{ℎ𝑘𝑖𝑗 (𝑥), −ℎ𝑘𝑖𝑗 (𝑥)}, 𝑘𝑖𝑗 =1 где 𝑓𝑘𝑖𝑗 (𝑥) = |(𝑋𝑖 − 𝑋𝑗 )𝑇 (𝑋𝑖 − 𝑋𝑗 ) − 𝑑2𝑖𝑗 |, ℎ𝑘𝑖𝑗 (𝑥) = (𝑋𝑖 − 𝑋𝑗 )𝑇 (𝑋𝑖 − 𝑋𝑗 ) − 𝑑2𝑖𝑗 . Здесь 𝑘𝑖𝑗 — порядковый номер пары (𝑖, 𝑗) такой, что 1 ≤ 𝑖 < 𝑗 ≤ 𝑚 и 𝑑𝑖𝑗 ̸= 0, а 𝑀 — количество таких пар. Выпишем в явном виде формулы для вычисления компонент градиентов функций ℎ𝑘𝑖𝑗 при 𝑘𝑖𝑗 ∈ 1 : 𝑀 𝜕ℎ𝑘𝑖𝑗 = 0𝑛 , 𝜕𝑋𝑝 ∀𝑝 ∈ 1 : 𝑚, 𝑝 ̸= 𝑖, 𝑝 ̸= 𝑗, 𝜕ℎ𝑘𝑖𝑗 = 2(𝑋𝑖 − 𝑋𝑗 ), 𝜕𝑋𝑖 𝜕ℎ𝑘𝑖𝑗 = −2(𝑋𝑖 − 𝑋𝑗 ). 𝜕𝑋𝑗 Аналогично предыдущему случаю, во всех градиентах ∇ℎ𝑘𝑖𝑗 отбросим первые 𝑛 компонент, соответствующие частным производным по компонентам точки 𝑋1 . Таким образом, ∇ℎ𝑘𝑖𝑗 будет иметь размерность 𝑁 = 𝑛(𝑚−1). 13

Теперь, зная градиенты гладких функций, не составит труда вычислить гиподифференциал функции 𝐹2 (𝑥) пользуясь формулами гиподифференциального исчисления (5)–(7). Важно отметить, что при увеличении числа атомов 𝑚 количество элементов в 𝑑𝐹2 (𝑥) будет расти экспоненциально. Это будет происходить в силу того, что максимальное число слагаемых в функционале 𝐹2 равняется числу ненулевых элементов 𝑑𝑖𝑗 при 1 ≤ 𝑖 < 𝑗 ≤ 𝑚, и оценивается сверху числом (𝑚−1)𝑚. Тогда согласно (6), чтобы сформировать гиподифференциал, потребуется порядка 𝑂(2(𝑚−1)𝑚 ) операций. Очевидно, что по возможности нужно избежать этой ситуации. Согласно алгоритму гиподифференциального спуска, гиподифференциал 𝑑𝐹2 (𝑥) понадобится лишь на этапе проверки необходимого условия минимума и поиска направления спуска, что сводится к нахождению минимального по норме гипоградиента. Однако, в силу устройства функции 𝐹2 , задачу минимизации (16) можно свести к задаче квадратичного программирования. Покажем это. Гиподифференциалы гладких функций ℎ𝑘𝑖𝑗 и −ℎ𝑘𝑖𝑗 (для краткости вместо 𝑘𝑖𝑗 далее будем писать просто 𝑘 ) при 𝑘 ∈ 1 : 𝑀 равны {︃(︃ 𝑑ℎ𝑘 (𝑥) = )︃}︃ 0 ∇ℎ𝑘 (𝑥) {︃(︃ , 𝑑(−ℎ𝑘 (𝑥)) = 0 )︃}︃ . −∇ℎ𝑘 (𝑥) Отсюда, согласно гиподифференциальному исчислению, получаем гиподифференциал функции 𝑓𝑘 (𝑥), 𝑘 ∈ 1 : 𝑀 (21) 𝑑𝑓𝑘 (𝑥) = co{𝐴𝑘 , 𝐵𝑘 }, где (︃ {︃ 𝐴𝑘 = 𝑑ℎ𝑘 (𝑥) + {︃ 𝐵𝑘 = −𝑑ℎ𝑘 (𝑥) + (︃ ℎ𝑘 (𝑥) − 𝑓𝑘 (𝑥) )︃}︃ 0𝑁 −ℎ𝑘 (𝑥) − 𝑓𝑘 (𝑥) 0𝑁 14 {︃(︃ = )︃}︃ ℎ𝑘 (𝑥) − 𝑓𝑘 (𝑥) ∇ℎ𝑘 (𝑥) {︃(︃ = )︃}︃ −ℎ𝑘 (𝑥) − 𝑓𝑘 (𝑥) −∇ℎ𝑘 (𝑥) , )︃}︃ .

Так как гиподифференциал 𝑑𝑓𝑘 (𝑥) это выпуклое множество, любую его точку 𝑣𝑘 ∈ 𝑑𝑓𝑘 (𝑥) можно представить в виде выпуклой комбинации 𝑣𝑘 = 𝛼𝑘 𝐴𝑘 + (1 − 𝛼𝑘 )𝐵𝑘 , где 𝛼𝑘 ∈ [0, 1]. С другой стороны гиподифференциал 𝑑𝐹2 (𝑥) = 𝑀 ∑︀ 𝑑𝑓𝑘 (𝑥) тоже вы- 𝑘=1 пуклое множество и любую его точку 𝑣 ∈ 𝑑𝐹2 (𝑥) можно представить как 𝑣= 𝑀 ∑︁ 𝑣𝑘 = 𝑘=1 𝑀 ∑︁ (︀ )︀ 𝛼𝑘 𝐴𝑘 + (1 − 𝛼𝑘 )𝐵𝑘 , 𝑘=1 где 𝛼𝑘 ∈ [0, 1]. Таким образом, задачу (16) свели к задаче квадратичного программирования 1 ‖𝑣‖2 → min . 2 𝑣∈𝑑𝐹2 (𝑥) Везде в программе, где необходимо решить задачу квадратичного программирования, использовали встроенную в MATLAB функцию quadprog. Задачу одномерной минимизации будем решать аналогично предыдущему случаю (см. пункт 3.1), используя функцию fmincon. В качестве правила останова используем каждый из критериев (18) (20) по очереди и сведем в единую таблицу количество итераций прошедших до выхода. Минимизация функционала 𝐹3 (𝑥) Гиподифференциал 𝑑𝐹3 (𝑥) легко вычисляется по формулам (5) и (7) и равен {︃ 𝑑𝐹3 (𝑥) = co 𝑑𝑓𝑘 (𝑥) + (︃ 𝑓𝑘 (𝑥) − 𝐹3 (𝑥) 0𝑁 )︃ ⃒ }︃ ⃒ ⃒ ⃒ 𝑘 ∈ 1:𝑀 , ⃒ (22) где 𝑑𝑓𝑘 (𝑥) находятся из (21). Заметим, что на этот раз число элементов в 𝑑𝐹3 (𝑥) не превышает (𝑚−1)𝑚. Перепишем формулу (22) в виде 𝑑𝐹3 (𝑥) = co{𝑍1 , . . . , 𝑍𝑠 }, 15 𝑠 = 2𝑀,

где 𝑍𝑖 ∈ R𝑁 +1 , при 𝑁 = 𝑛(𝑚−1), и находятся по формулам (︃ 𝑍2𝑘−1 = ℎ𝑘 (𝑥) − 𝐹3 (𝑥) ∇ℎ𝑘 (𝑥) )︃ (︃ , 𝑍2𝑘 = −ℎ𝑘 (𝑥) − 𝐹3 (𝑥) −∇ℎ𝑘 (𝑥) )︃ , 𝑘 ∈ 1 : 𝑀. Согласно алгоритму гиподифференциального спуска, после того, как найден гиподифференциал 𝑑𝐹3 (𝑥) необходимо решить задачу минимизации (16). На этот раз для ее решения эффективно воспользоваться МДМметодом, алгоритм которого описан в главе 2. Заметим, что МДМ-метод плохо сходится, когда начало координат принадлежит 𝑑𝐹3 (𝑥). В связи с этим в [5] рекомендуется делать следующее. Вначале решим задачу линейного программирования в пространстве векторов 𝑊 = (𝛼1 , . . . , 𝛼𝑠 , 𝑢): найдем min 𝑢 при условиях 𝑠 ∑︁ − 𝑖=1 𝑠 ∑︁ (𝑖) 𝑘 ∈ 1 : 𝑁 +1, (𝑖) 𝑘 ∈ 1 : 𝑁 +1, 𝛼𝑖 𝑧𝑘 − 𝑢 ≤ 0, 𝛼𝑖 𝑧𝑘 − 𝑢 ≤ 0, 𝑖=1 𝑠 ∑︁ 𝛼𝑖 = 1; 𝛼𝑖 ≥0, 𝑖 ∈ 1 : 𝑠; 𝑢 ≥ 0, 𝑖=1 (𝑖) (𝑖) где 𝑍𝑖 = (𝑧1 , . . . , 𝑧𝑁 +1 ) — точки из 𝑑𝐹3 (𝑥). (0) (0) Обозначим через 𝑊0 = (𝛼1 , . . . , 𝛼𝑠 , 𝑢(0) ) решение этой задачи. Если 𝑢(0) = 0, то начало координат содержится в гиподифференциале, и алгоритм гиподифференциального спуска завершен. Если же 𝑢(0) > 0, то для нахождения ближайшей к началу координат точки гиподифференциала, можно воспользоваться МДМ-методом, взяв в качестве начального приближения точку 𝑍0 = 𝑠 ∑︁ (0) 𝛼𝑖 𝑍𝑖 , 𝑖=1 (0) (0) где вектор (𝛼1 , . . . , 𝛼𝑠 ) составлен из первых 𝑠 компонент вектора 𝑊0 . 16

На практике условие останова (20) заменим эквивалентным условием |𝑢(0) | < 𝜀, (23) а задачу линейного программирования решим, воспользовавшись функцией linprog. Задачу одномерной минимизации решаем аналогично предыдущим случаям (см. пункт 3.1). В качестве правила останова, помимо критерия (23), будем использовать критерии (18) и (19), и для каждого критерия сведем в единую таблицу количество итераций прошедших до выхода. 17

4 Численные эксперименты Все рассмотренные выше алгоритмы реализованы на языке MATLAB. Для каждого из функционалов 𝐹1 , 𝐹2 , 𝐹3 были написаны сценарии, реализующие соответствующие алгоритмы минимизации. На вход каждой функции передаются четыре параметра: матрица расстояний 𝐷, начальная точка 𝑥0 , точность вычислений 𝜀 и максимальное количество итераций 𝑁 𝑢𝑚𝐼𝑡𝑒𝑟. На выходе получаем точку минимума соответствующего функционала. Рассмотрим несколько примеров с разными входными параметрами. Во всех примерах положим 𝑁 𝑢𝑚𝐼𝑡𝑒𝑟 = 200 и 𝜀 = 10−4 . В критериях (12), (13), (18) и (19) будем использовать непосредственно 𝜀 = 10−4 , а в критериях (14), (20) и (23) используем 𝜀 = 10−2 (в связи со спецификой критериев). Заметим, что произвольно сформированная матрица 𝐷 не подойдет для нашей задачи, т.к. в общем случае по ней нельзя восстановить корректную конфигурацию молекулы. Будем формировать матрицу 𝐷 на основе координат атомов реально существующей молекулы. Выберем произвольную молекулу из Protein Data Bank (банк данных 3-D структур белков и нуклеиновых кислот). Пусть это будет молекула с идентификатором 1PTQ, состоящая из 404 атомов. Во избежание больших размерностей в примерах будем использовать не все атомы молекулы, а лишь некоторые из них. По известным координатам этих атомов сформируем полную матрицу расстояний 𝐷. В качестве начальной точки 𝑥0 выберем координаты (𝑥𝑖 + 𝛿𝑥 , 𝑦𝑖 + 𝛿𝑦 , 𝑧𝑖 + 𝛿𝑧 ), где 𝑖 ∈ 1 : 𝑚, (𝑥𝑖 , 𝑦𝑖 , 𝑧𝑖 ) — известные координаты атомов, а 𝛿𝑥 , 𝛿𝑦 , 𝛿𝑧 — случайные величины, подчиняющиеся стандартному нормальному распределению. Тем самым мы смоделировали ситуацию, когда известна полная матрица расстояний и некоторая начальная точка близкая к 𝑥* — оптимальной точке, доставляющей минимум всем трем функционалам. 18

Пример 4.1. Случай полной матрицы, 𝑚 = 10. Результат работы алгоритмов минимизации функционалов 𝐹1 , 𝐹2 и 𝐹3 приведен в таблицах 1, 2 и 3 соответственно. 𝑘 𝑓 (𝑥𝑘+1 ) |𝑓 (𝑥𝑘 ) − 𝑓 (𝑥𝑘+1 )| ‖𝑥𝑘 − 𝑥𝑘+1 ‖ ‖∇𝑓 (𝑥𝑘 )‖ 1 8260.2707511790 8260.2707511790 13.2544287096 6040.3248898236 5 453.9828111224 198.4693130093 0.6888006248 603.8107615152 20 47.5342117678 4.9514857255 0.0774235571 127.1388025029 50 2.2129246930 0.2605837340 0.0172937807 30.1074781694 100 0.0411799823 0.0028448104 0.0018859290 3.0164233547 150 0.0015528047 0.0000992959 0.0003559610 0.5590147346 200 0.0001296652 0.0000057771 0.0000953830 0.1213709544 Таблица 1: минимизация 𝐹1 в случае полной матрицы, 𝑚 = 10 𝑘 𝑓 (𝑥𝑘+1 ) |𝑓 (𝑥𝑘 ) − 𝑓 (𝑥𝑘+1 )| ‖𝑥𝑘 − 𝑥𝑘+1 ‖ ‖¯ 𝑧𝑘 ‖ 1 144.2262555770 144.2262555770 2.5913798097 127.6774733645 5 3.8337048298 15.2395353420 0.8233148493 12.3558759441 10 0.0001334643 0.0007932857 0.0000277705 0.0009267500 17 0.0000100027 0.0000097699 0.0000003087 0.0000002328 Таблица 2: минимизация 𝐹2 в случае полной матрицы, 𝑘 𝑓 (𝑥𝑘+1 ) |𝑓 (𝑥𝑘 ) − 𝑓 (𝑥𝑘+1 )| 1 12.4168327175 12.4168327175 5 2.1041535953 10 13 ‖𝑥𝑘 − 𝑥𝑘+1 ‖ 𝑚 = 10 ‖𝑢(𝑘) ‖ 2.9441990362; 4.8797177542 0.7673412804 0.5271780222 0.4223796526 1.0855159106 0.3827402898 0.9943154442 0.1381279400 0.0064073252 0.0000000105 0.0000002778 0.0063155860 Таблица 3: минимизация 𝐹3 в случае полной матрицы, 𝑚 = 10 Видим, что метод градиентного спуска довольно медленно сходился и превысил максимальное число итераций. Метод гиподифференцильного спуска, напротив, сошелся быстро: за 17 и 13 итераций для 𝐹2 и 𝐹3 соответственно. Максимальной точности по значению функции в точке минимума достиг функционал 𝐹2 . Время работы градиентного спуска составило 33.66с., гиподифференциального — 5.93с. и 8.60с. при минимизации 𝐹1 и 𝐹2 соответственно. 19

По всем параметрам наилучшим методом оказался метод гиподифференциального спуска минимизирующий функционал 𝐹1 . Пример 4.2. Случай неполной матрицы, 𝑚 = 6. Обнулив некоторые случайно выбранные компоненты корректной матрицы, получим разреженную матрицу ⎛ 0 1, 464 0 0 2, 458 0 ⎞ ⎟ ⎜ ⎟ ⎜1, 464 0 1, 525 2, 411 0 2, 591 ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ 0 1, 525 0 1, 232 2, 459 0 ⎟. 𝐷=⎜ ⎟ ⎜ 2, 411 1, 232 0 0 4, 620⎟ ⎜ 0 ⎟ ⎜ ⎜2, 458 0 2, 459 0 0 1, 519⎟ ⎠ ⎝ 0 2, 591 0 4, 620 1, 519 0 Результат работы алгоритмов приведен в таблицах 4, 5 и 6. Аналогично предыдущему примеру медленней всех отработал метод градиентного спуска — за 13.8с. Время работы алгоритмов гиподифференциального спуска при минимизации 𝐹1 и 𝐹2 составило 1.57с. и 2.55с. соответственно. В случае разреженной матрицы расстояний алгоритм оптимизирующий функционал 𝐹2 показал лучшие результаты. Алгоритм градиентного спуска — худшие. 𝑘 𝑓 (𝑥𝑘+1 ) |𝑓 (𝑥𝑘 ) − 𝑓 (𝑥𝑘+1 )| ‖𝑥𝑘 − 𝑥𝑘+1 ‖ ‖∇𝑓 (𝑥𝑘 )‖ 1 0.3295250778 0.3295250778 0.1699607237 72.6172899711 20 0.0168285291 0.0018310064 0.0036965003 0.9901334297 50 0.0007277416 0.0000776693 0.0007665760 0.2031338062 100 0.0000133201 0.0000010312 0.0000824726 0.0258778964 111 0.0000069185 0.0000003078 0.0000943951 0.0096069820 Таблица 4: минимизация 𝐹1 в случае неполной матрицы, 𝑚=6 𝑘 𝑓 (𝑥𝑘+1 ) |𝑓 (𝑥𝑘 ) − 𝑓 (𝑥𝑘+1 )| ‖𝑥𝑘 − 𝑥𝑘+1 ‖ ‖¯ 𝑧𝑘 ‖ 1 0.2397530326 0.2397530326 0.2831620837 4.5378906335 4 0.0000007361 0.0000005537 0.0000003076 0.0000012898 Таблица 5: минимизация 𝐹2 в случае неполной матрицы, 20 𝑚=6

𝑘 𝑓 (𝑥𝑘+1 ) |𝑓 (𝑥𝑘 ) − 𝑓 (𝑥𝑘+1 )| ‖𝑥𝑘 − 𝑥𝑘+1 ‖ ‖𝑢(𝑘) ‖ 1 0.0714320502 0.0714320502 0.2924567705 1.4933117648 4 0.0014927833 0.0000003999 0.0000023039 0.0014862652 Таблица 6: минимизация 𝐹3 в случае неполной матрицы, 𝑚=6 Ориентируясь на приведенные примеры, можно сделать вывод, что метод гиподифференциального спуска, в сравнении с методом градиентного спуска, сходится за гораздо меньшее количество итераций, работает быстрее, а также, в случае функционала 𝐹2 , показывает наилучшую точность. 21

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

ЛИТЕРАТУРА [1] Dong Q., Wu Z. A linear time algorithm for solving the molecular distance geometry problem with exact inter-atomic distances. // Journal of Global Optimization 22: pp. 365-375, 2002 [2] Grosso A., Locatelli M., Schoen F. Solving molecular distanse geometry problem by global optimization algorithms // Springer Science, 2007. pp. 23-37 [3] Thi L., An H. Solving large scale molecular distance geometry problems by a smoothing technique via the gaussian transform and D.C. programming // Journal of global optimization 27: 375-397, 2003. [4] Глебов Н. И., Кочетов Ю. А., Плясунов А. В. Методы оптимизации. // Новосиб. ун-т. Новосибирск, 2000. 105с. [5] Демьянов В. Ф., Малозёмов В. Н. Введение в минимакс. // Издательство «Наука», 1972. 368с. [6] Демьянов В. Ф., Васильев Л. В. Недифференцируемая оптимизация. // М.: Наука, 1981. 384с. [7] Демьянов В. Ф., Рубинов А. М. Основы негладкого анализа и квазидифференциальное исчисление. // М.: Наука, 1990. 432c. [8] Малозёмов В. Н. МДМ-методу — 40 лет. // Семинар «DHA& CAGD». Избранные доклады. 10 декабря 2011. 23

Приложение Основной файл сценарий. 1 close all ; clc ; 2 getpdb ( ' 1 ptq ' , ' T o F i l e ' , ' 1 p t q _ f i l e . p d b ' ) ; % скачиваем молекулу 1PTQ 3 p d b s t r u c t = pdbread ( ' 1 p t q _ f i l e . p d b ' ) ; % считываем файл в структуру 4 atom = pdbstruct.Model.Atom ; % извлекаем поле атом 5 S = 1 0 ; % число считываемых атомов 6 XYZ = z e r o s ( S , 3 ) ; % массив координат 7 f o r i =1:S 8 XYZ( i , 1 ) = atom ( i ) .X ; 9 XYZ( i , 2 ) = atom ( i ) .Y ; 10 XYZ( i , 3 ) = atom ( i ) . Z ; 11 end 12 [m, n ] = s i z e (XYZ) ; 13 x0 = z e r o s (m, n ) ; 14 s h i f t = XYZ( 1 , : ) ; % сдвиг координат в ноль 15 f o r i =1:m x0 ( i , : ) = XYZ( i , : ) - s h i f t ; 16 17 18 19 end D = z e r o s (m, m) ; f o r i =1:m- 1 f o r j=i +1:m 20 D( i , j ) = s q r t ( ( x0 ( i , : ) - x0 ( j , : ) ) * ( x0 ( i , : ) - x0 ( j , : ) ) ' ) ; 21 end 22 23 24 25 26 end y0 = x0 ; % всем координатам добавляем случайную погрешность f o r i =2: l e n g t h ( y0 ) 27 y0 ( i , 1 ) = y0 ( i , 1 )+ randn ; 28 y0 ( i , 2 ) = y0 ( i , 2 )+ randn ; 29 y0 ( i , 3 ) = y0 ( i , 3 )+ randn ; 30 end 31 e p s = 0 . 0 0 0 1 ; i t e r = 2 0 0 ; % точность и 32 x1 = GradientDescentF1 ( D, y0 , eps , i t e r ) ; 33 D1 = CheckAnswer ( x1 , D) ; % градиентный спуск для F1 число итераций 34 t i t l e ( 'Молекула найденная минимизацией F1' ) ; 35 figure 36 x2 = H y p o d i f f e r e n t i a l D e s c e n t F 2 ( D, y0 , eps , i t e r ) ; 37 D2 = CheckAnswer ( x2 , D) ; % гиподифференциальный спуск для F2 38 t i t l e ( 'Молекула найденная минимизацией F2' ) ; 39 figure ; 40 x3 = H y p o d i f f e r e n t i a l D e s c e n t F 3 ( D, y0 , eps , i t e r ) ; 41 D3 = CheckAnswer ( x3 , D) ; % гиподифференциальный спуск для F3 42 t i t l e ( 'Молекула найденная минимизацией F3' ) ; 24

Минимизация 𝐹1 (𝑥) методом градиентного спуска. 1 function y = GradientDescentF1 ( D, x0 , eps1 , IterNum ) 2 % Минимизация функционала F1 ( x ) методом градиентного с п у с к а . 3 % D - матрица расстояний , x0 - начальная точка , 4 % e p s 1 - точность , IterNum - число итераций 5 [m, n ] = s i z e ( x0 ) ; 6 d i a = max(max(D) ) * 1 0 ; % 10 - ти кратный " диаметр " молекулы 7 xk = x0 ; i t e r = 0 ; 8 % флаги критериев останова 1 , 2 , 3 9 flag_g = 0 ; flag_y = 0 ; flag_x = 0 ; 10 f = z e r o s ( 2 , 1 ) ; % значение функции на смежных итерациях 11 w h i l e i t e r < IterNum 12 iter = iter + 1; 13 gk = GradF1 ( xk , D ) ; 14 fhnd_Func1 = @( t ) Func1 (D, xk , gk , t ) ; % h a n d l e функция 15 o p t s = o p t i m s e t ( ' D i s p l a y ' , ' none ' ) ; 16 [ tk , f (mod( i t e r , 2 ) +1) ] = fmincon ( fhnd_Func1 , d i a / 2 , [ ] , [ ] , [ ] , [ ] , 0 , ... % задача одномерной минимизации dia , [ ] , o p t s ) ; d i a = tk * 1 0 ; prev_xk = xk ; % предыдущая точка 17 18 20 xk ( 2 :m, : ) = xk ( 2 :m, : ) - gk * tk ; % следующее приближение f k = f (mod( i t e r , 2 ) +1) ; 21 n f k = abs ( f ( 1 ) - f ( 2 ) ) ; 22 i f n f k < e p s 1 && f l a g _ y 19 ̸= 1 d i s p ( [ ' 1 ) Ответ по норме приращения функции. Итерация № ' ... 23 num2str ( i t e r ) ' . F1 ( x * ) = ' num2str ( f k ) ] ) ; flag_y = 1 ; 24 25 end 26 ngk = Norma ( gk ) ; % норма градиента 27 i f ngk < e p s 1 * 100 && f l a g _ g ̸= 1 d i s p ( [ ' 3 ) Ответ по норме г р а д и е н т а . Итерация № ' ... 28 num2str ( i t e r ) ' . F1 ( x * ) = ' num2str ( f k ) ] ) ; flag_g = 1 ; 29 30 end 31 Δ 32 nxk = Norma ( Δ ) ; 33 if = prev_xk - xk ; nxk < e p s 1 && f l a g _ x ̸= 1 d i s p ( [ ' 2 ) Ответ по норме приращения аргумента. Итерация № ' ... 34 num2str ( i t e r ) ' . F1 ( x * ) = ' num2str ( f k ) ] ) ; flag_x = 1 ; 35 36 end 37 i f f l a g _ x * f l a g _ y * f l a g _ g == 1 break ; % если в с е критерии сработали - выход 38 end 39 40 end 25

41 if iter IterNum ≥ d i s p ( [ ' Алгоритм превысил ' , num2str ( IterNum ) , ' итераций. ... 42 F1 ( x * ) = ' num2str ( min ( f ) ) ] ) ; 43 end 44 y = xk ; 45 end Минимизация 𝐹2 (𝑥) методом гиподифференциального спуска. 1 function y = H y p o d i f f e r e n t i a l D e s c e n t F 2 ( D, x0 , eps1 , IterNum ) 2 % Минимизация функционала F2 ( x ) методом гиподифференциального с п у с к а . 3 % D - матрица расстояний , x0 - начальная точка , 4 % e p s 1 - точность , IterNum - число итераций 5 [m, n ] = s i z e ( x0 ) ; % m число атомов , n размерность 6 N = n * (m - 1 ) ;%число неизвестных 7 d i a = max(max(D) ) * 2 ; 8 M = nnz (D) ; %количество ненулевых элементов в D 9 xk = x0 ; i t e r = 0 ; 10 % флаги критериев останова 1 , 2 , 3 11 flag_g = 0 ; flag_y = 0 ; flag_x = 0 ; 12 f = zeros (2 ,1) ; 13 w h i l e i t e r < IterNum i t e r = i t e r +1; k=1; 14 15 h i = z e r o s (M, 1 ) ; 16 hi_grad = z e r o s (M, N) ; 17 f o r i =1:m- 1 % вычисляем градиенты f o r j=i +1:m 18 if 19 D( i , j ) ̸= 0 20 h i ( k ) = h i j ( xk ( i , : ) , xk ( j , : ) , D( i , j ) ) ; 21 hi_grad ( k , : ) = Gradhi ( xk ( i , : ) , xk ( j , : ) , n * ( i - 1 ) - ( n - 1 ) , ... n * ( j - 1 ) - ( n - 1 ) , N) ; 22 k=k+1; end 23 end 24 25 end 26 h = max( hi , - h i ) ; % формируем гиподифференциал 27 A = [ hi - h , hi_grad ] ; B = [ - hi - h , - hi_grad ] ; 28 % находим ближайшее расстояние до начала координат 29 g i p o g r a d i e n t = DistanceToZero (A, B) ' ; 30 ngk = norm ( g i p o g r a d i e n t ) ; % норма гипоградиента 31 z = g i p o g r a d i e n t ( 2 :N+1) ; 32 g = - z /norm ( z ) ; % направление спуска 33 gk = z e r o s (m- 1 , n ) ; 34 k=1; 35 f o r i =1:m- 1 36 f o r j =1:n 26

gk ( i , j ) = g ( k ) ; 37 k=k+1; 38 end 39 40 end 41 fhnd_Func2 = @( t ) Func2 (D, xk , gk , t , M) ; % h a n d l e функция 42 % одномерная минимизация 43 o p t s = o p t i m s e t ( ' D i s p l a y ' , ' none ' ) ; 44 [ tk , f (mod( i t e r , 2 ) +1) ] = fmincon ( fhnd_Func2 , d i a / 2 , [ ] , [ ] , [ ] , [ ] , 0 , ... dia , [ ] , o p t s ) ; 45 d i a = tk * 5 ; 46 fk = 47 prev_xk = xk ; % предыдущее приближение 48 xk ( 2 :m, : ) = xk ( 2 :m, : ) + gk * tk ; % следующее приближение nyk = abs ( f ( 1 ) - f ( 2 ) ) ; 49 f (mod( i t e r , 2 ) +1) ; i f nyk < e p s 1 && f l a g _ y 50 ̸= 1 d i s p ( [ ' 1 ) Ответ по норме приращения функции. Итерация № ' ... 51 num2str ( i t e r ) ' . F2 ( x * ) = ' num2str ( f k ) ] ) ; flag_y = 1 ; 52 53 end 54 i f ngk < e p s 1 * 100 && f l a g _ g ̸= 1 d i s p ( [ ' 3 ) Ответ по норме гипоградиента. Итерация № ' ... 55 num2str ( i t e r ) ' . F2 ( x * ) = ' num2str ( f k ) ] ) ; flag_g = 1 ; 56 57 end 58 nxk = Norma ( prev_xk - xk ) ; 59 if nxk < e p s 1 && f l a g _ x ̸= 1 60 d i s p ( [ ' 2 ) Ответ по норме приращения аргумента. Итерация № ' ... 61 flag_x = 1 ; num2str ( i t e r ) ' . F2 ( x * ) = ' num2str ( f k ) ] ) ; 62 end 63 i f f l a g _ x * f l a g _ y * f l a g _ g == 1 break ; % если в с е критерии сработали - выход 64 end 65 66 end 67 if iter ≥ IterNum d i s p ( [ ' Алгоритм превысил ' , num2str ( i t e r ) , ' итераций. F2 ( x * ) = ' ... 68 num2str ( f k ) ] ) ; 69 end 70 fclose ( fid ) ; 71 y = xk ; 72 end 27

Минимизация 𝐹3 (𝑥) методом гиподифференциального спуска. 1 function y = H y p o d i f f e r e n t i a l D e s c e n t F 3 ( D, x0 , eps1 , IterNum ) 2 % Минимизация функционала F3 ( x ) методом гиподифференциального с п у с к а . 3 % D - матрица расстояний , x0 - начальная точка , 4 % e p s 1 - точность , IterNum - число итераций 5 [m, n ] = s i z e ( x0 ) ; % m число атомов , n размерность 6 M = nnz (D) ; % количество ненулевых элементов в D 7 N = 8 xk = x0 ; i t e r = 0 ; 9 % флаги критериев останова 1 , 2 , 3 n * (m - 1 ) ; d i a = max(max(D) ) ; 10 flag_g = 0 ; flag_y = 0 ; flag_x = 0 ; 11 f = zeros (2 ,1) ; 12 w h i l e i t e r < IterNum i t e r = i t e r + 1 ; k=1; 13 14 15 16 h i = z e r o s (M, 1 ) ; hi_grad = z e r o s (M, N) ; % вычисляем градиенты гладких функций f o r i =1:m- 1 f o r j=i +1:m 17 if 18 D( i , j ) ̸= 0 19 h i ( k ) = h i j ( xk ( i , : ) , xk ( j , : ) , D( i , j ) ) ; 20 hi_grad ( k , : ) = Gradhi ( xk ( i , : ) , xk ( j , : ) , n * ( i - 1 ) - ( n - 1 ) , ... n * ( j - 1 ) - ( n - 1 ) , N) ; 21 k=k+1; end 22 end 23 24 end 25 F3 = max(max( hi , - h i ) ) ; % формируем гиподифференциал 26 A = [ hi - F3 , hi_grad ] ; B = [ - hi - F3 , - hi_grad ] ; 27 Gipo = [A; B] ' ; 28 % проверяем лежит ли ноль в гиподифференциале 29 30 [ F l a g N u l l I n s i d e , z0 , u ] = Z e r o I n s i d e ( Gipo , 10 * e p s 1 ) ; i f F l a g N u l l I n s i d e == 1 31 d i s p ( [ ' 3 ) Ответ по норме z0 найден на 32 if iter < 2 ' num2str ( i t e r ) ' итерации. ' ] ) ; f ( 1 ) = Func2 (D, xk , xk ( 2 :m, : ) , 0 , M) ; 33 34 end 35 d i s p ( [ 'F3 ( x * ) = ' num2str ( min ( f ) ) ] ) ; 36 end 37 g i p o g r a d i e n t = MDM_method( Gipo , z0 ) ; % МДМ- метод 38 ngk = norm ( g i p o g r a d i e n t ) ; % норма гипоградиента 39 z = g i p o g r a d i e n t ( 2 :N+1) ; 40 g = - z /norm ( z ) ; % направление спуска 41 gk = z e r o s (m- 1 , n ) ; 42 k=1; 43 f o r i =1:m- 1 28

f o r j =1:n 44 45 gk ( i , j ) = g ( k ) ; 46 k=k+1; end 47 48 end 49 % задача одномерной оптимизации 50 fhnd_Func3 = @( t ) Func3 (D, xk , gk , t , M) ; % h a n d l e функция 51 o p t s = o p t i m s e t ( ' D i s p l a y ' , ' none ' ) ; 52 [ tk , f (mod( i t e r , 2 ) +1) ] = ... fmincon ( fhnd_Func3 , d i a / 2 , [ ] , [ ] , [ ] , [ ] , 0 , dia , [ ] , o p t s ) ; 53 d i a = tk * 5 ; 54 fk = 55 prev_xk = xk ; % точка на предыдущей итерации 56 xk ( 2 :m, : ) = xk ( 2 :m, : ) + gk * tk ; % следующее приближение nyk = abs ( f ( 1 ) - f ( 2 ) ) ; 57 f (mod( i t e r , 2 ) +1) ; i f nyk < e p s 1 && f l a g _ y 58 ̸= 1 d i s p ( [ ' 1 ) Ответ по норме приращения функции. Итерация № ' ... 59 num2str ( i t e r ) ' . F3 ( x * ) = ' num2str ( f k ) ] ) ; flag_y = 1 ; 60 61 end 62 i f ngk < e p s 1 * 100 && f l a g _ g ̸= 1 d i s p ( [ ' 3 ) Ответ по норме гипоградиента. Итерация № ' ... 63 num2str ( i t e r ) ' . F3 ( x * ) = ' num2str ( f k ) ] ) ; flag_g = 1 ; 64 65 end 66 Δ 67 nxk = Norma ( Δ ) ; 68 if = prev_xk - xk ; nxk < e p s 1 && f l a g _ x ̸= 1 d i s p ( [ ' 2 ) Ответ по норме приращения аргумента. Итерация № ' ... 69 num2str ( i t e r ) ' . F3 ( x * ) = ' num2str ( f k ) ] ) ; flag_x = 1 ; 70 71 end 72 i f f l a g _ x * f l a g _ y * f l a g _ g == 1 break ; % если в с е критерии сработали - выход 73 end 74 75 76 end 77 if iter ≥ IterNum d i s p ( [ ' Алгоритм превысил ' , num2str ( i t e r ) , ' итераций. F3 ( x * ) = ' ... 78 num2str ( f k ) ] ) ; 79 end 80 y = xk ; 81 end 29

Функция 𝐹1 (𝑥). 1 f u n c t i o n y = Func1 ( D, x0 , g , t ) 2 % функицонал F1 ( x ) . D - матрица расстояний , 3 % x0 - начальная точка , g - направление спуска , 4 % t - параметр одномерной минимизации 5 [m, n ] = s i z e ( x0 ) ; 6 x = z e r o s (m, n ) ; 7 x ( 2 :m, : ) = x0 ( 2 :m, : ) - t * g ; 8 f i = z e r o s (m- 1 , 1 ) ; 9 k=1; 10 f o r i =1:m- 1 f o r j=i +1:m 11 i f D( i , j ) ̸= 0 12 f i ( k ) = ( ( x ( i , : ) - x ( j , : ) ) * ( x ( i , : ) - x ( j , : ) ) ' - D( i , j ) ^2) ^ 2 ; k=k+1; 13 14 end 15 end 16 17 end 18 y = sum ( f i ) ; 19 end Функция 𝐹2 (𝑥). 1 f u n c t i o n y = Func2 ( D, x0 , g , t , M) 2 % функицонал F2 ( x ) . D - матрица расстояний , 3 % x0 - начальная точка , g - направление спуска , 4 % t - параметр одномерной минимизации , 5 % M - число ненулевых элементов в матрице D 6 [m, n ] = s i z e ( x0 ) ; 7 x = z e r o s (m, n ) ; 8 x ( 2 :m, : ) = x0 ( 2 :m, : ) + t * g ; 9 f i = z e r o s (M, 1 ) ; 10 k=1; 11 f o r i =1:m- 1 f o r j=i +1:m 12 i f D( i , j ) ̸= 0 13 f i ( k ) = ( ( x ( i , : ) - x ( j , : ) ) * ( x ( i , : ) - x ( j , : ) ) ' - D( i , j ) ^2) ; k = k+1; 14 15 end 16 end 17 18 end 19 y = sum (max( f i , - f i ) ) ; 20 end 30

Функция 𝐹3 (𝑥). 1 f u n c t i o n y = Func3 ( D, x0 , g , t , M) 2 % функицонал F3 ( x ) . D - матрица расстояний , 3 % x0 - начальная точка , g - направление спуска , 4 % t - параметр одномерной минимизации , 5 % M - число ненулевых элементов в матрице D 6 [m, n ] = s i z e ( x0 ) ; 7 x = z e r o s (m, n ) ; 8 x ( 2 :m, : ) = x0 ( 2 :m, : ) + t * g ; 9 f i = z e r o s (M, 1 ) ; 10 k=1; 11 f o r i =1:m- 1 f o r j=i +1:m 12 i f D( i , j ) ̸= 0 13 14 f i ( k ) = ( ( x ( i , : ) - x ( j , : ) ) * ( x ( i , : ) - x ( j , : ) ) ' - D( i , j ) ^2) ; 15 k = k+1; end 16 end 17 18 end 19 y = max(max( f i , - f i ) ) ; 20 end Градиент функции 𝐹1 (𝑥). 1 f u n c t i o n y = GradF1 ( x0 , D ) 2 % функция вычисляющая градиент функции F1 ( x ) 3 % x0 - точка , D - матрица расстояний 4 [m, n ] = s i z e ( x0 ) ; 5 Phi = z e r o s ( n*m, m- 1 ) ; % матрица частных производных функций 6 f o r i =1:m- 1 f o r j=i +1:m 7 i f D( i , j ) ̸= 0 8 buf = ( ( x0 ( i , : ) - x0 ( j , : ) ) * ( x0 ( i , : ) - x0 ( j , : ) ) ' - D( i , j ) ^2) * ( ... 9 x0 ( i , : ) - x0 ( j , : ) ) ' ; 10 % формируем диагональные элементы матрицы 11 12 Phi ( n* i - ( n - 1 ) : n* i , i ) = Phi ( n* i - ( n - 1 ) : n* i , i ) + 4 * buf ; % формируем элементы под диагональю 13 Phi ( n* j - ( n - 1 ) : n* j , i ) = -4 * buf ; end 14 end 15 16 end 17 y = z e r o s (m- 1 , n ) ; 18 % первые n строк нас не интересуют , т . к . они отвечают частным ... производным по 19 % компонентам 1 - го атома , который мы зафиксировали в начале координат 31

20 k=n+1; 21 f o r i =2:m f o r j =1:n 22 23 y ( i - 1 , j ) = sum ( Phi ( k , : ) ) ; 24 k=k+1; end 25 26 end 27 end Функция ℎ𝑘𝑖𝑗 (𝑥). 1 f u n c t i o n y = h i j ( xi , xj , d ) 2 % функция , вычисляющая значение функции h i j в точках xi , x j . 3 % d - ( i , j ) -ый элемент матрицы D 4 y = ( xi - x j ) * ( xi - x j ) ' - d ^ 2 ; 5 end Градиент функции ℎ𝑘𝑖𝑗 (𝑥). 1 f u n c t i o n g = Gradhi ( xi , xj , i , j , N ) 2 % функция , вычисляющая градиенты гладких функций h i ( xi , x j ) 3 % xi , x j - аргументы hi , i , j - индексы аргументов 4 % N - общее число неизвестных задачи 5 6 n = l e n g t h ( x i ) ; % размерность 7 g = z e r o s (N, 1 ) ; 8 i f i <1 % В силу специфики функций h i j - ых , на i - ую позицию ставим значения : g ( j : j+n - 1 , 1 ) = 2 * xj ' ; 9 10 else 11 g ( i : i+n - 1 , 1 ) = 2 * ( x i - x j ) ' ; 12 % на j - ую позицию ставим : 13 g ( j : j+n - 1 , 1 ) = 2 * ( x j - x i ) ' ; 14 end 15 % Остальные компоненты градиента g равны нулю 16 end 32

Поиск минимального расстояния от начала координат до выпуклого множества. 1 f u n c t i o n vmin = DistanceToZero ( A, B ) 2 % функция поиска минимального расстояния от гиподифференциала 3 % до начала координат методом квадратичного программирования. [ A, B ] - ... гиподифференциал 5 % задаем матрицу H для задачи квадратичного программирования (xT*H* x ) AB = A-B ; H = AB*AB' ; 6 % формируем остальные параметры функции quadprog 7 Bs = sum (B) ; f = Bs *AB' ; 8 [M, N1 ] = s i z e (A) ; 4 9 % верхние и нижние ограничения 10 l b = z e r o s (M, 1 ) ; ub = o n e s (M, 1 ) ; 11 o p t s = o p t i m s e t ( ' L a r g e S c a l e ' , ' on ' , ' D i s p l a y ' , ' none ' ) ; 12 % рещаем задачу квадратичного программирования 13 [ alpha ] = quadprog (H, f , [ ] , [ ] , [ ] , [ ] , lb , ub , [ ] , o p t s ) ; 14 vmin = Bs + alpha '*AB; 15 end Проверка на принадлежность начала координат выпуклому множеству. 1 f u n c t i o n [ N u l l I n s i d e , z0 , u ] = Z e r o I n s i d e ( A1 , e p s 2 ) 2 % функция , проверяющая принадлежность нуля выпуклому множеству A1 3 % если ноль принадлежит , то возвращает N u l l I n s i d e =1, z0=0 . 4 % иначе N u l l I n s i d e =0, z0 - начальное приближение для МДМ метода. 5 % u - искомый параметр метода 6 [ size_m , s i z e _ n ] = s i z e (A1) ; 7 % формируем последний столбец коэффициентов вектора u 8 last_column = - o n e s ( size_m * 2 , 1 ) ; 9 b = z e r o s ( size_m * 2 , 1 ) ; % столбец b 10 % матрица A для ограничений неравенств 11 A = [ A1 ; -A1 ] ; 12 A = [ A, last_column ] ; 13 % вектор - столбец коэффициентов целевой функции ( минимизируем только 14 % последнюю компоненту u ) 15 f = z e r o s ( s i z e _ n +1, 1 ) ; 16 f ( s i z e _ n +1, 1 ) = 1 ; 17 % матрица Aeq для ограничений равенств 18 Aeq = o n e s ( 1 , s i z e _ n +1) ; 19 Aeq ( 1 , s i z e _ n +1) = 0 ; 20 beq = 1 ; 21 % ограничения снизу l b 22 ≤ x ≤ ub l b = z e r o s ( s i z e _ n +1, 1 ) ; 33

% линейное программирование. используем симплекс метод 23 24 o p t i o n s = o p t i m s e t ( ' Simplex ' , ' on ' , ' D i s p l a y ' , ' none ' ) ; 25 [ alfa_u ] = l i n p r o g ( f , A, b , Aeq , beq , lb , [ ] , [ ] , o p t i o n s ) ; 26 z0 = z e r o s ( size_m , 1 ) ; 27 u = abs ( alfa_u ( s i z e _ n +1) ) ; 28 i f u < eps2 NullInside = 1; 29 else 30 31 NullInside = 0; 32 f o r i =1: s i z e _ n z0 = z0 + alfa_u ( i ) *A1 ( : , i ) ; 33 end 34 end 35 36 end МДМ-метод. 1 f u n c t i o n [ v_min ] = MDM( A, v0 ) 2 % МДМ метод. Поиск ближайшей к нулю точки симплекса 3 % A - симплекс , v0 - начальное приближение 4 e p s 2 = 0 . 0 0 1 ; % точность МДМ- метода 5 vk = v0 ; 6 [ N, M] = s i z e (A) ; 7 % вектор из единиц для составления линейной системы с ограничениями 8 a = o n e s ( 1 , M) ; 9 % решаем систему , учитывая , что p_1 + p_2 + . . . + p_m = 1 & p_i 10 11 12 13 ≥ 0 o p t s = o p t i m s e t ( 'TolX' , 0 . 0 0 0 0 0 1 ) ; p = l s q n o n n e g ( [ A; a ] , [ vk ; 1 ] , [ ] , o p t s ) ; N = 1 0 0 0 ; % максимальное число итераций f o r i t e r =1:N 14 M_plus= f i n d ( p ) ; % формируем носитель вектора p 15 % ищем индекс i ' , на котором реалиуется max( a_i , v_k) , по всем i из ... M_plus 16 max_ai_vk = vk '*A( : , M_plus ) ; 17 min_ai_vk = vk '*A; % находим максимум и запоминаем индекс i 1 18 19 [ max_i1 , i1_with_hat ] = max( max_ai_vk ) ; 20 i 1 = M_plus ( i1_with_hat ) ; 21 22 % находим минимум и запоминаем индекс i 2 [ min_i2 , i 2 ] = min ( min_ai_vk ) ; 23 % вычисляем дельта и проверяем условие останова 24 Δ 25 = max_i1 - min_i2 ; if Δ < eps2 26 v_min = vk ; 27 break ; 28 end 34

29 Δ _a1a2 = A( : , i 1 ) -A( : , i 2 ) ; 30 v_hat = vk - p ( i 1 ) *Δ _a1a2 ; 31 norma_a1a2 = 32 tk_hat = 33 i f tk_hat ≥ 1 34 tk = 1 ; Δ Δ _a1a2 '*Δ _a1a2 ; / ( p ( i 1 ) * norma_a1a2 ) ; else 35 tk = tk_hat ; 36 37 end 38 vk = vk + tk * ( v_hat - vk ) ; 39 % вычисляем p 40 pi1 = p( i1 ) ; pi2 = p( i2 ) ; 41 f o r i =1: l e n g t h ( p ) if i 42 ̸= i 1 && i ̸= i 2 p ( i )=p ( i ) ; 43 else 44 i f i == i 1 45 p ( i ) = ( 1 - tk ) * p i 1 ; 46 47 end 48 i f i == i 2 p ( i ) = p i 2 + tk * p i 1 ; 49 end 50 end 51 end 52 53 end 54 v_min = vk ; 55 % следующее приближение end Вычисление нормы Евклида. 1 2 f u n c t i o n a = Norma ( x ) % разворачиваем матрицу x ( [ m, n ] ) в вектор и вычисляем е г о норму 3 [m, n ] = s i z e ( x ) ; 4 a = 0; 5 f o r i =1:m f o r j =1:n 6 a = a + x ( i , j ) ^2; 7 end 8 9 end 10 a = sqrt (a) ; 11 end 35

Восстановление матрицы 𝐷 по заданному вектору 𝑥. 1 f u n c t i o n D_ans = CheckAnswer ( x , D ) 2 % функция , восстановливающая матрицу D по заданному вектору x 3 % и рисующая молекулу. 4 [m, n ] = s i z e ( x ) ; 5 D_ans = z e r o s (m, m) ; 6 k=1; 7 f o r i =1:m- 1 f o r j=i +1:m 8 D_ans ( i , j ) = s q r t ( ( x ( i , : ) - x ( j , : ) ) * ( x ( i , : ) - x ( j , : ) ) ' ) ; 9 i f D( i , j ) ̸= 0 10 11 bonds ( k , : ) = [ i , j ] ; 12 k=k+1; end 13 end 14 15 end 16 % рисуем ответ 17 i f n == 2 % двумерный случай p l o t ( x ( : , 1 ) , x ( : , 2 ) , 'o' ) ; 18 f o r i =1:m 19 text (x( i ,1) , x( i ,2) , [' 20 ' num2str ( i ) ] , ' c o l o r ' , ' r e d ' ) ; end 21 f o r k=1: l e n g t h ( bonds ) 22 l i n e ( [ x ( bonds ( k , 1 ) , 1 ) x ( bonds ( k , 2 ) , 1 ) ] , [ x ( bonds ( k , 1 ) , 2 ) , ... 23 x ( bonds ( k , 2 ) , 2 ) ] ) ; end 24 g r i d on ; 25 26 end 27 i f n==3 % трехмерный случай p l o t 3 ( x ( : , 1 ) , x ( : , 2 ) , x ( : , 3 ) , 'o' ) ; 28 f o r i =1:m 29 text (x( i ,1) , x( i ,2) , x( i ,3) , [' 30 ' num2str ( i ) ] , ' c o l o r ' , ' r e d ' ) ; end 31 f o r k=1: l e n g t h ( bonds ) 32 l i n e ( [ x ( bonds ( k , 1 ) , 1 ) , x ( bonds ( k , 2 ) , 1 ) ] , [ x ( bonds ( k , 1 ) , 2 ) , ... 33 x ( bonds ( k , 2 ) , 2 ) ] , [ x ( bonds ( k , 1 ) , 3 ) , x ( bonds ( k , 2 ) , 3 ) ] ) ; 34 end 35 g r i d on ; 36 end 37 end 36

Рецензии:

Авторизуйтесь, чтобы добавить рецензию

- у работы пока нет рецензий -

Отзывы:

Авторизуйтесь, чтобы оставить отзыв