Построение единой программной среды для решения задач глобальной оптимизации

МИНИСТЕРСТВО ОБЩЕГО И ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ РФ

АЛТАЙСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
Математический факультет
Кафедра геоинформационных технологий

Построение единой программной среды для решения задач глобальной оптимизации

Выполнил студент

5 курса, 441группы

Ахмеров Р.Р.

_________________________

(подпись)

Научный руководитель

ст. преподаватель

Жилин С.И.

_________________________

(подпись)

Допустить к защите Дипломная работа защищена

Зав. кафедрой "___"______________ 1999 г.

Поляков Ю.А., д.т.н, доцент ОЦЕНКА ________________

_______________________ Председатель ГАК

(подпись) _________________________

Ф.И.О.

"____"__________1999г. _________________________

(подпись)

Барнаул 1999
Реферат

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

Работа содержит 6 рисунков, две таблицы и одно приложение. Общий объем работы составляет 39 страниц.

Содержание
Содержание 2
Введение 3

Обзор основных методов глобальной оптимизации 4

Стохастические методы 5

Моделируемый Отжиг 5

Методы кластеризации 6

Детерминированные методы 8

Метод ветвей и границ 8

Функции ограничения 10

Ограничение, основанное на константах Липшица 11

Интервальная арифметика
11

Обобщенная интервальная арифметика
12

Аффинная арифметика
13

Смешанная интервально-аффинная арифметика 16

Управление ошибкой округления 17

Оценка границ значений функций посредством

интервального разложения в ряд Тейлора
18

Схемы разбиения
19

Ускорение интервальных методов ветвей и границ 20

Back-Boxing 20

Аффинное оценивание элементарных функций 21

Аффинная оценка функции 23

Аффинная оценка функции
e

x
24

Представление аффинных форм на ЭВМ 25

Общие подформулы 25

Практическая реализация алгоритмов
27

Автоматический анализ функций 28

Эффективное распределение памяти 30

Структурирование результатов 32

Численные результаты 32
Заключение 36
Литература 37
Приложение 38
Введение

Общая задача, рассматриваемая в данной работе, состоит в поиске минимума ¦* некоторой целевой функции ¦:W®R, определенной на компактном множестве WÍRn
и множества всех минимизирующих переменных W*(¦) = {x
* Î W: ¦(x
*) = ¦*}. Поскольку в общем случае компактное множество может иметь сложную структуру здесь полагается, что W является параллелепипедом [x
1
, 1
]´[ x
2
, 2
]´…´ [ x
n
, n
] Í Rn
.

Глобальная оптимизация без ограничений – одна из математических задач, рассматривается до сих пор и является очень сложной. Множество исследователей занимались решением этой задачи и некоторыми из них были приведены случаи NP
-полных задач глобальной оптимизации. Точное доказательство NP-
полноты очень трудоемко, но сложность решения не требует доказательства. В качестве примеров таких нелинейных задач, в которых поиск минимума (минимума или максимума) необходим для практики, можно привести задачу поиска направления и продолжительности радиационного облучения раковых опухолей, поиска состояния наименьшей энергии протеиновой молекулы, и задачу определения правильного соотношения составляющих смеси для оптимального протекания химической реакции.

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

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

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

Обзор основных методов глобальной оптимизации

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

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

· Значения, которые аппроксимируют глобальный минимум с математически рассчитанной максимальной ошибкой.

· Функция не ограничена.

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

Стохастические методы

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

Алгоритм 1
: Случайный поиск

Инициализация
:
¦min
= ¥, xmin
= any_number, и i
= 0.

шаг 1
:
xi
= случайное число.

шаг 2
:
если ¦( xi
) < ¦min
установить xmin
= xi
, ¦min
=¦( xi
).

шаг
3:
увеличить i
на 1 и перейти к шагу 1.

Если ¦ хотя бы кусочно непрерывна, то при i
® ¥, ¦min
® min ¦ на выбранной области. Этот процесс будет сходиться к глобальному минимуму не только на Rn
, но и на любом множестве машинно-представимых чисел в случае реализации этой процедуры на компьютере.

Этот метод был изучен несколькими исследователями. Они показали, что если функция вычисляется в точках, равномерно распределенных в области W, то поиск наименьшего значения функции таким методом сходится к глобальному минимальному значению ¦*
с вероятностью 1.

Моделируемый Отжиг

Моделируемый отжиг (simulated annealing - SA) - это другой стохастический метод, дающий хорошие результаты в целочисленной оптимизации, и, с недавнего времени, в общей глобальной оптимизации.

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

Алгоритм 2
: Простой Моделируемый Отжиг

Инициализация
:
Положить x0
= 0, и i
= 0.

Шаг 1
:
Выбрать y
в некоторой окрестности xi
.

Шаг 2
:
Если ¦(y
) < ¦( xi
) установить xi+1
= y
, перейти к шагу 4.

Шаг 3
:
Если (равномерно распределенное случайное число)<
, установить xi+1
= y
.

Шаг 4:
Увеличить i
на 1 и перейти к шагу 1.

здесь c(
i)
есть температура в момент времени i
, и ¦ минимизируемая функция. В общем случае ожидается, что c(
i)
убывает до некоторого c
¥
¹ 0 при i
® ¥.

Для дискретного случая разработка этого алгоритма достаточно ясна. Достаточно только определить окрестность xi
и функцию остывания c(
i)
. Когда SA применяется к непрерывной задаче, идея окрестности даже более не ясна. Для непрерывного SA можно вначале определить направление поиска. Это делается путем случайного выбора точки на единичной гиперсфере с центром в xi
, которая и дает нам нужное направление. После этого алгоритм выбирает случайную длину шага в этом направлении. Очевидным недостатком этой процедуры является то, что она не использует какую-либо локальную информацию (такую как дифференцируемость). Для учета этой проблемы предлагается случайным образом сочетать следующих два метода. Либо xi
выбирается случайно, либо xi
выбирается путем применения нескольких шагов поиска локального минимума для xi-1
.

Методы кластеризации

Методы кластеризации (clustering methods) относятся к стохастическим методам. Метод кластеризации стартует с некоторой точки, случайно выбранной из области W, и затем создает группы или кластеры ближайших точек, соответствующих основной области притяжения. Кластер C
есть множество точек, которое соответствуют области притяжения и имеет следующие свойства:

· C
содержит ровно один локальный минимум x*
.

· Любой алгоритм спуска, стартующий из точки x
Î C,
сходится к точке x*
.

· Первые два свойства подразумевают, что C
может содержать только одну стационарную точку.

Например, рассмотрим функцию (x
2
– 1)2
, которая имеет глобальный минимум в точках x
= ±1. Для нее существует два кластера. Первый кластер это {x
| -< x
< 0}, и второй кластер это {x
| 0< x
< } (см. рис. 1). Любой метод спуска, стартующий в кластере 1, гарантированно сходится к точке x
= -1. Аналогично, метод спуска, стартующий в кластере 2, гарантированно сходится к точке x
= 1.

Cluster 1 Cluster 2 Рис. 1. Кластеры функции (x 2 -1)2

-1.5 -1 0 1 1.5

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

После проведения кластеризации из всех полученных локальных минимумов выбирается точка с наименьшим значением функции.

Детерминированные методы

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

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

Метод ветвей и границ

Метод ветвей и границ (branch and bound) принадлежит классу детерминированных методов и может применяться для решения широкого круга задач, в которых информация о минимизируемых функциях мала или отсутствует совсем. Он может быть легко использован для задач с вогнутыми или выпуклыми областями, вогнутыми или выпуклыми целевыми функциями и с функциями, для которых недоступна информация об их непрерывности. В тех случаях, когда для работы алгоритма требуется какая-либо специфическая информация об исследуемой функции, существует две возможности: либо пользователь предоставляет эту информацию на входе, либо алгоритм сам проводит аналитическое исследование функции, вычисляя нужные производные, определяя области непрерывности и т.д. Очевидно, что на практике более удобно работать с вычислителем, основанном на последнем подходе, но бывают случаи, когда получение информации (например, вычисление константы Липшица) аналитическим путем является задачей, плохо поддающейся, или совсем не поддающейся алгоритмизации. В таких случаях получение всей необходимой информации возлагается на пользователя.

Метод ветвей и границ, подобно моделируемому отжигу, первоначально успешно применялся для решения задач целочисленной оптимизации. После этого он был адаптирован для непрерывного случая и для решения задачи поиска min¦(x
) при x
ÎWÍRn
принял следующий канонический вид:

Алгоритм 3
: Метод ветвей и границ

Инициализация
:
Установить i
= 1, z*
= + infinite
, и список = W.

Шаг 1
:
Извлечь область Wi
из списка.

Шаг 2
:
Если Wi
достаточно мала, занести ее в результирующий список и перейти к шагу 1.

Шаг 3
:
Вычислить z
-
, z
+
такие что z
-
£ ¦(x
)|W
i
£ z
+
.

Шаг 4
:
Если z
-
> z*
удалить область Wi
и перейти к шагу 1.

Шаг 5
:
Если z
+
£ z*
установить z*
= z
+
.

Шаг 6
:
Разбить Wi
на k
равных частей и добавить их в список.

Шаг 7
:
Перейти на шаг 1.

Шаг 8
:
Вернуть список результирующих областей.

Алгоритм 3 разбивает исходную область W на все более и более мелкие части, отбрасывая те, наличие глобального минимума в которых невозможно.

Этот простой алгоритм имеет несколько недостатков. Во-первых, границы значений функции на Wi
, вычисленные на шаге 3 потенциально не точны. Если оценка функции проведена с плохой точностью, то алгоритм может потратить большое количество времени, пытаясь отбросить области, не содержащие глобального минимума, на которых целевая функция является достаточно плоской. Во-вторых, на шаге 6, Wi
разбивается на k
равных частей, причем не принимается в расчет какая-либо информация об локально/глобальном поведении функции на полученных частях. Более интеллектуальный алгоритм должен принять эту информацию во внимание, для определения того, где нужно резать область. Например, возьмем прямоугольник [-1, 1]´[-1, 1] и ¦(x
) = x
2
. Этот прямоугольник можно разбить на четыре части, каждая с углом в (0, 0). Это разбиение не очень хорошее, поскольку каждый прямоугольник все еще содержит глобальный минимум функции. В то же время, если прямоугольник разрезать на части

[-e
, e
]´ [-e
, e
], [e
, 1]´ [-1, 1], [-1, -e
]´ [-1, 1], [-e
, e
]´ [e
, 1], [-e
, e
]´ [-1, -e
],

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

В-третьих, подобласти могут отбрасываться только на основе информации о значениях функции. Если же функция хотя бы один раз непрерывно дифференцируема и показано, что градиент этой функции не может принимать нулевое значение на Wi
, целиком лежащей (т.е. не имеющей общих граничных точек) в исходном W, то Wi
не содержит даже локального минимума. Аналогично, подобласть может быть отброшена, если матрица Гессе (гессиан) ¦ ни в одной точке этой подобласти не является положительно определенной.

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

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

Функции ограничения

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

Ограничение, основанное на константах Липшица

Потенциально недорогая, но неточная оценка функции может быть дана следующим выражением

BL
,
W
¦ = ¦(x
*
) – L
w
(W) £ ¦(x
) £ ¦(x
*
) + L
w
(W) = BL
,
W
¦,

где x*
Î W,
а L
есть константа Липшица для функции ¦ на W, и w
(W) служит верхней границей диаметра области W. Такая оценка имеет то преимущество, что она проста для вычисления, когда известна константа Липшица. Дополнительно она имеет следующее свойство

lim (BL
,
W
¦ - BL
,
W
¦) = 0 при w
(W)® 0.

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

Интервальная арифметика

Интервальная арифметика (IA
– interval arithmetic
) – естественная тех­ника для оценки границ изменения функции. IA
была введена Муром [2] с целью улучшения надежности численных вычислений. С тех пор она успешно приме­нялась для решения многих вычислительных задач, таких как обыкновенные дифференциальные уравнения, системы линейных уравнений и глобальная оптимизация.

В IA
вещественное число x
представимо парой чисел с плавающей точ­кой (a
, b
), соответствующих интервалу X
= [a
, b
], гарантированно содержащему x
, т.е. такому, что a
£ x
£ b
. Таким образом, IA
обеспечивает не только оценку для значения x
, но и говорит о том, насколько хороша эта оценка. Реальная сила IA
состоит в том, что мы можем оперировать интервалами, как если бы они были числами, и получать надежные оценки для результатов числовых вычислений.

Простые формулы легко получены для выполнения примитивных ариф­метических операций над интервалами. Интервальные расширения для сложных функций могут быть вычислены путем составления этих примитивных формул в том же порядке, в каком они составлены при вычислении непосредственно самой функции. Другими словами, любой алгоритм для вычисления функции, использующий примитивные операции может быть легко (и автоматически) пе­реведен в алгоритм для вычисления интервального расширения для этой же функции. Это особенно удобно выполнять в языках программирования, кото­рые поддерживают перегрузку операторов, но может быть выполнено в любом языке программирования, либо вручную, либо при помощи предварительной компиляции. Поскольку также относительно легко найти интервальное расши­рение для элементарных трансцендентных функций, таких как sin
, cos
, log
и exp
, класс функций, для которых интервальное расширение может быть легко (и ав­томатически) вычислено гораздо больше, чем класс рациональных полиноми­альных функций. Эти наблюдения иногда объединяют в Фундаментальную Теорему Интервальной Арифметики: каждая вычислимая функция имеет интервальное расширение, т.е. для любой вещественной функции ¦, заданной алгоритмом, существует интервальная функция F
, такая, что F
(X
) Ê ¦(C
) = {¦(x
): x
ÎX
} для каждого параллелепипеда X
в области определения ¦. Причем lim F
(X
) = ¦(x
), при сужении параллелепипеда X
к любой точки x
из области определения ¦.

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

Данный недостаток IA
особенно сильно проявляется в длинных вычислительных цепях, где часто наблюдается “взрыв ошибки”: поскольку оценка продвигается ниже по цепи, относительная точность вычисленных интервалов уменьшается экспоненциально, и они вскоре становятся слишком широкими для использования. Подходы к решению проблемы зависимости в IA
включают использование центрированных форм, обобщенной интервальной арифметики Хансена, и более новой – аффинной арифметики [5].

Обобщенная интервальная арифметика

Работая над проблемой зависимости в IA
, Хансен [11] изобрел новую вычислительную модель, которая называется обобщенной интервальной арифметикой (GIA
– generalized interval arithmetic
). В этой модели, каждая входная переменная xi
представлена интервалом, как в IA
, и каждая промежуточная величина z
, вычисленная в процессе оценки ¦(x1,…,
xn
) представлена как линейная комбинация входных переменных:

= z0
+ z1
x1
+…+ zn
xn
,

где каждый коэффициент zi
сам является интервалом. Только коэффициенты zi
записываются для каждой величины z
; переменные xi
подразумеваются. Т.о. GIA
представляет величины многочленами первой степени с интервальными коэффициентами.

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

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

= [5, 7] x1
+ [3, 4] x2
,

= [4, 6] x1
+ [4, 5] x2
,

где переменные x1
и x2
имеют диапазоны [-10, 10]. Из этих данных мы можем вывести, что y
и z
содержатся в интервале [-110, 110]. В GIA
модели разность
=
–
оценивается

= [-1, 3] x1
+ [-2, 0] x2
.

Это означает, что w
находится в [-30, 30], в то время как IA
дает [-220, 220].

Аффинная арифметика

Другая модель, применяемая к проблеме зависимости – аффинная арифметика (AA
– affine arithmetic
) [4]. В AA
каждая величина x
представлена аффинной формой

= x0
+ x1
e1
+ x2
e2
+ … + xn
en
,

где x0
– центр формы, xi
(i
= 1
,..,n
) – известные вещественные коэффициенты (записанные как числа с плавающей точкой), называемые частичными отклонениями, и ei
– символьные переменные, называемые символами шума, чьи значения могут произвольно изменяться в интервале [-1, 1]. Т.о. AA
поверхностно подобна GIA
в том, что обе они представляют величины полиномами первой степени от символьных переменных. Тем не менее, как мы увидим ниже, число членов используемых в AA
растет при выполнении вычисления, в то время как число членов используемых в GIA
фиксировано и равно числу входных переменных.

Аффинные формы неявно представляют частичные зависимости между операндами: когда две аффинные формы совместно используют общие символы шума, величины, ими представленные, по крайней мере, частично зависят одна от другой. Как и в GIA
, принятие таких зависимостей во внимание позволяет AA
обеспечивать намного более точные оценки диапазонов, чем в IA
, особенно в длинных вычислительных цепях. Другая основная польза, которая связана с этим, – то, что ошибка аффинного представления квадратично зависит от размера W, а не линейно.

AA
используется в диапазонном анализе следующим образом. Вначале преобразовываем все входные интервалы к аффинным формам. Затем оперируем этими аффинными формами с помощью AA
для вычисления желаемой функции. Окончательно, преобразовываем результат назад в интервал. Преобразования между IA
и AA
представлениями очевидны. Данный интервал
= [a
, b
] представляет некоторую величину x
эквивалентной аффинной формой =
x0
+ xk
ek
,
где x0
= (b+a)
/2, xk
= (b-a)
/2. Поскольку входные интервалы обычно представляют независимые переменные, то для каждого входного интервала должен быть использован новый символ шума ek
. Обратно, значение величины, представленной аффинной формой
= x0
+ x1
e1
+ x2
e2
+ … + xn
en
, гарантированно будет лежать в интервале [x0
-
x
, x0
+
x
], где x
= ||
|| = å|xi
|. Заметим, что [x0
-
x
, x0
+
x
] – самый маленький интервал, содержащий все возможные значения
, в предположении, что каждый ei
изменяется независимо в интервале [-1, 1].

Для оценки функции ¦ в AA
, мы должны заменить каждую примитивную операцию Ф
, которая появляется в выражении ¦ на эквивалентную операцию над аффинными формами, как это сделано в IA
для интервалов. Для линейных операций z
= Ф
(x, y)
соответствующая аффинная форма
может быть выражена точно, как линейная комбинация символов шума, присутствующих в формах
и
. Т.е. если

= x0
+ x1
e1
+ x2
e2
+ … + xn
en
,

= y0
+ y1
e1
+ y2
e2
+ … + yn
en
,
и a
ÎR
, то:

±

= (x0

±
y0
) + (x1
±
y1
)
e1
+ (x2

±
y2
)
e2
+ … + (xn

±
yn
)
en
,

a
= (
a
x0
) + (
a
x1
)
e1
+ (
a
x2
)
e2
+ … + (
a
xn
)
en
,

±
a
= (x0
+
a
) + x1
e1
+ x2
e2
+ … + xn
en

Заметим, что разность
–
равна 0, – полезное свойство, которое не верно в IA
(или даже в GIA
). Отсутствие таких тривиальных сокращений в IA
является одним из источников взрыва ошибки.

Когда примитивная операция Ф
нелинейна, значение
не может быть выражено непосредственно как линейная комбинация ei
. В этом случае, мы подбираем хорошую линейную аппроксимацию для Ф
(“хорошую” в чебышевском смысле минимизации максимальной ошибки), и добавляем дополнительный член zk
ek
для представления ошибки этой аппроксимации:

= z0
+ z1
e1
+ z2
e2
+ … + zn
en
+ zk
ek
.

Здесь ek
– новый шумовой символ и zk
– верхняя граница ошибки аппроксимации. Выбор линейной аппроксимации должен принимать во внимание диапазоны операндов, которые подразумеваются аффинными формами
и
, и, возможно, также их взаимосвязь. Отметим, что в отличие от GIA
, число членов, используемых в аффинных формах, растет во время оценки формулы. (При необходимости, лишние члены могут быть объединены и заменены новыми шумовыми символами ценою некоторой потери точности)

Важно заметить, что для дифференцируемой операции Ф
граница ошибки z
k
хорошей чебышевской аппроксимации зависит квадратично от диапазонов операндов ||
|| и ||
||. Т.о., если функция ¦ использует только дифференцируемые операции, то различие между действительным диапазоном ¦ и диапазоном, вычисленном в AA
стремится к нулю в относительном смысле при сужении области. Напротив, внутренние ошибки IA
(следствие неучтенных взаимосвязей) зависят линейно от размеров области; поэтому различие будет уменьшаться только в абсолютном смысле.

Адекватные аппроксимационные алгоритмы были разработаны для основных арифметических операций и трансцендентных функций [5]. Например, при умножении двух аффинных форм
и
мы можем использовать аппроксимацию

z0
= x0
y0
, zi
= x0
yi
+ y0
xi
, zk
= ||
||×||
||.

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

Как описывалось применительно к IA
, формулы для оценки примитивных функций в AA
могут быть автоматически объединены в формулы для произвольно сложных функций. Этот факт и возможность преобразования позволяют AA
просто заменить IA
.

В качестве простого примера того, как AA
обрабатывает проблему зависимости, рассмотрим оценку z = x(10 – x)
в AA
для x
в интервале [4, 6]:

= 5 + 1e1

10 –
= 5 - 1e1

=
(10 –
) = 25 + 0e1
+ 1e2

[
] = [25 - 1, 25 + 1] = [24,26].

Т.о. AA
вычисляет оценку очень близкую к [24, 25] – действительному диапазону z
. С другой стороны, IA
оценка есть [16, 36]. Если выражение x
(10 – x
) раскрыто в 10x – x
2
, то IA
оценка становится значительно хуже, [4, 44], в то время как оценка в AA
точна – [24, 25].

Смешанная интервально-аффинная арифметика

Преимущество аффинной арифметики над интервальной арифметикой всецело зависит от ее способности определить и использовать зависимости между операндами. Когда такая зависимость отсутствует, или случайно согласуется с консервативными предположениями IA
, AA
оценки могут быть шире, чем в IA
– даже не смотря на то, что аффинная арифметика более точна, чем интервальная.

Объяснением этого парадоксального факта служит то, что IA
пытается окружить график каждой элементарной функции ортогональным параллелепипедом минимального вертикального размера, в то время как AA
ищет наклонный параллелепипед минимального объема. Например, рассмотрим операцию возведения в квадрат z = x
2
, для x
Î[1, 3]. IA
дает (узкий) диапазон [1, 9] для z
; а AA
представляет x
как
=
2+1e1
, и производит
=
4.5 + 4e1
+ 0.5e2
, чей диапазон [0, 9]. Заметим, что IA
говорит только то, что пара (x
, z
) лежит в прямоугольнике площади 2 * 8 = 16, в то время как результат AA
ограничен параллелограммом площади 2 * 1 = 2.

Этот недостаток оценок диапазонов в AA,
к сожалению, имеет тенденцию проявляться, когда в функции ¦ преобладают члены, подобные квадрату, оцениваемые возле своего минимума – именно там, где нам нужны точные нижние границы.

Эта проблема может быть решена путем объединения IA
с AA
. В этой смешанной модели каждая величина x
представлена и аффинной формой
и стандартным интервалом [x
]. Каждая операция z = Ф
(x, y)
выполняется в обеих вычислительных моделях . Интервальная часть [z
] результата – это множество, являющееся пересечением диапазонов, вычисленных в IA
и в AA
; и AA-
процедура использует интервалы [x
] и [y
] при выборе чебышевской аппроксимации для Ф
. Таким образом, IA
и AA
находятся во взаимодействии на каждом шаге: явные интервалы, полученные IA,
приводят к лучшим выборам для самых лучших линейных аппроксимаций и к меньшим значениям для члена ошибки z
k
; в то время как аффинные формы в AA
позволяют сокращать связанные ошибки, и, следовательно, приводить к более узким интервалам. В частности, смешанная AAIA
(affine arithmetic – interval arithmetic
) схема вычисляет неотрицательные диапазоны для Ф
(x) = x
2
, без потери информации о зависимости, обеспеченной AA
.

Управление ошибкой округления

Когда используется точная арифметика, IA
и AA
всегда вычисляют математически допустимые границы для диапазона значений последовательности операций. Эти границы могут включать реальный диапазон, но они всегда правильны. Однако IA
и AA
, будут выполняться на цифровом компьютере, который использует арифметику с плавающей точкой, которая подвержена ошибкам округления. Но все же IA
и AA
способны произвести математически допустимые границы на любом вычислительном устройстве, поддерживающем контроль над округлением. Наиболее современные вычислительные архитектуры обеспечивают контроль над ошибками округления, реализуя IEEE стандарт чисел с плавающей точкой [4].

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

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

В IA
много дополнительной работы требуется для направленного округления переменных. Для AA
количество требуемой работы более чем удвоено для разных примитивных операций. К счастью, это не верно для алгоритма в целом, исключая случай, когда окончательная точность решения недостижима из-за машинной точности.

Оценка границ значений функций посредством

интервального разложения в ряд Тейлора

Любую n
+1 раз непрерывно дифференцируемую функцию ¦ в окрестности точки x
*
можно разложить в ряд Тейлора

¦(x
) = + Rn
(x
, x
),

где

Rn
(x
, x
) = ,

и x
лежит на отрезке, соединяющим точки x
и x
*
. Отсюда следует, что для любого x
из интервала X
значение ¦(x
) будет принадлежать интервалу

Tn
(X)
= +,

где вещественные операции заменены на соответствующие им интервальные операции и x
*
ÎX
. Tn
(X)
называется интервальным многочленом Тейлора n
-го порядка. Как и для других методов, для Tn
(X)
справедливо условие

Более того, если функция ¦ бесконечное число раз непрерывно дифференцируема, то интервальное разложение в ряд Тейлора имеет то свойство, что

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

Схемы разбиения

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

Способы определения того, как разбивать данный параллелепипед сильно варьируются. Параллелепипед можно делить пополам вдоль самого длинного ребра, делить пополам вдоль наиболее перспективного направления, или разбивать в соответствии с длиной градиента. Можно также разбивать на k
частей, где k
некоторое положительное число большее, либо равное 2. Выбор конкретного разбиения определяется из условия задачи, причем для теоретической сходимости алгоритма нужно гарантировать то, что длина при таком разбиении самого длинного ребра стремится к 0.

Ускорение интервальных методов ветвей и границ

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

Если целевая функция ¦ дифференцируема достаточное количество раз, то простой способ ускорить метод ветвей и границ состоит в использовании одного из многих быстрых методов локальной оптимизации [10]. Локальная оптимизация может быть совмещена с интервальными проверками для градиента или гессиана: если интервальная оценка для градиента Ñ
¦ на подобласти D не содержит 0, то D не может содержать локальный минимизатор. Если интервальная оценка для гессиана показывает, что он положительно определен на D, то D также не может содержать локальный минимизатор.

Back-Boxing

Back-Boxing
[1] – методика ускорения для ал

горитмов ветвей и границ, объединяющая процедуры быстрого локального поиска с результатом проверяющих методов для устранения и уменьшения большого объема выполненной работы, когда она выполняется вблизи локального минимума. Цель Back-Boxing
’а состоит в том, чтобы определить большие области в W, где целевая функция ¦ гарантированно имеет один локальный минимум. Если только такие области известны, их локальный минимум может быть вычислен эффективно и затем проверен, является ли он глобальным минимумом.

Back-Boxing
первоначально локализует минимум x
*
в некоторой подобласти D. Затем ищется наибольший параллелепипед B
ÍD, такой, что x
*
- уникальный локальный минимум для B
; область DB
затем помещается в список L для дальнейшей обработки. Поскольку мы теперь знаем, что параллелепипед B
имеет уникальный минимум, нам в дальнейшем не нужно делить B
, в попытке найти меньшее значение функции.

Аффинное оценивание элементарных функций

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

Как уже было сказано, чтобы оценить формулу в AA
мы должны заменить каждую из элементарных операций z
= ¦(x
, y
) на вещественных числах эквивалентной операцией =(,) над аффинными формами [4,7], где - процедура, вычисляющая аффинную форму z
= ¦(x
, y
), которая соответствует , .

По определению

x
= x
0
+ x
1
e1
+ … + x
n
en
(1)

y
= y
0
+ y
1
e1
+ …+ y
n
en
(2)

для некоторых (неизвестных) значений (e1
, …, en)
ÎUn
, где Un
есть декартово произведение n
вещественных интервалов [-1, 1]. Поэтому величина z
есть функция ei
, а именно

z
= ¦(x
, y
) = ¦(x
0
+ x
1
e1
+ … + x
n
en
, y
0
+ y
1
e1
+ …+ y
n
en
) = ¦*(e1
,…,en
) (3)

Процедура теперь должна заменить ¦*(e1
,…,en
) аффинной формой

= z
0
+ z
1
e1
+ …+ z
n
en,

которая сохраняет столько информации, сколько возможно из соотношений между x
, y
и z
, вытекающих из (1-3), но без использования других ограничений, которые не могут быть выведены из начальных данных.

В общем случае, когда ¦ является нелинейной операцией, функция ¦*(e1
,…,en
) = ¦(x
, y
) из (3) не может быть выражена как линейная комбинация ei
. В этом случае мы можем подобрать некоторую линейную функцию

¦a
(e1
,…,en
) = z
0
+ z
1
e1
+ …+ z
n
en
, (4)

которая аппроксимирует ¦*(e1
,…,en
) равномерно на области Un
, и затем добавить к ней дополнительный член z
k
ek
, представляющий ошибку аппроксимации. Таким образом,

z = ¦a
(e1
,…,en
) + z
k
ek
= z
0
+ z
1
e1
+ …+ z
n
en
+ z
k
ek.

Здесь ek
должен быть новым символом шума (отличным от всех других символов шума в одном вычислении) и величина z
k
должна быть верхней границей модуля разницы между ¦a
и ¦* для всех возможных значений e1
,…,en
; т.е.

max { | ¦*(e1
,…,en
) - ¦a
(e1
,…,en
) | : (e1
,…,en
)ÎUn
}

Заметим, что замена ¦* на ¦a
+ z
k
ek
частично отходит от основной цели AA
: символ шума ek
принимается независимым от e1
,…,en
, в то время как фактически это функция от них. Любые последующие операции, использующие z
в качестве входного значения, не имеют сведений о взаимосвязи между ek
и e1
,…,en
и поэтому могут возвращать аффинную форму, которая менее точна, чем необходимо.

Для минимизации потерянной информации мы должны выбрать коэффициенты z
0
, z
1
,…, z
n
таким образом, чтобы сделать новый член ошибки как можно меньше.

Другими словами ¦a
должна быть многочленом первой степени, который наилучшим образом аппроксимирует ¦* на Un
в смысле Чебышева минимизации максимальной ошибки.

Отметим одно свойство, которое окажется нам полезным в дальнейшем. Пусть нам надо найти наилучшую линейную аппроксимацию для функции ¦(x
(e1
,…,en
)), где x
(e1
,…,en
) = x
0
+ x
1
e1
+ … + x
n
en
, тогда ¦*(e1
,…,en
) = ¦(x
0
+ x
1
e1
+ … + x
n
en
) аппроксимируется на гиперкубе Un
некоторой функцией ¦a
(e1
,…,en
) = z
0
+ z
1
e1
+ …+ z
n
en
. Заметим, что функция ¦*(e1
,…,en
) – константа на любой гиперплоскости из Un
, ортогональной вектору (x
1
,…, x
n
). Нетрудно показать, что наилучшая (по Чебышеву) линейная аппроксимация ¦* также должна обладать этим свойством. Т.е. нам нужно рассматривать только аппроксимации ¦a
(e1
,…,en
) вида

¦a
(e1
,…,en
) = a + b = a( x
0
+ x
1
e1
+ … + x
n
en
) + b.

Также легко показать, что значения a и b, которые минимизируют максимальную ошибку ¦a
– ¦* есть в точности коэффициенты наилучшей линейной аппроксимации функции ¦(t
) функцией at +
b на интервале []. Таким образом, задача аппроксимации исходной функции n
переменных редуцируется в задачу аппроксимации функции одной переменной. Можно показать (см. [16]), что для функции одной переменной прямая at +
b будет наилучшей линейной аппроксимацией на отрезке [a
, b
], если

a = , b = (5)

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

Аффинная оценка функции

Итак, рассмотрим аффинную оценку функции [8]. Нам нужно аппроксимировать функцию ¦*(e1
,…,en
) = на гиперкубе Un
некоторой функцией ¦a
(e1
,…,en
) = z
0
+ z
1
e1
+ …+ z
n
en
и затем добавить к последней дополнительный член z
k
ek
для учета ошибки аппроксимации, где ek
вновь созданный символ шума.

Как показано ранее, нам нужно рассматривать только аппроксимации вида ¦a
(e1
,…,en
) = a + b = a( x
0
+ x
1
e1
+ … + x
n
en
) + b.

Если [a
, b
] – интервал [], то оптимальные чебышевские коэффициенты a и b есть

a =

b =

и максимальная ошибка аппроксимации d =

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

После того как a, b и d получены, мы можем вычислить

z
0
= ax
0
+ b

z
i
= ax
i
(i = 1,…,n)

z
k
= d

Этот анализ подразумевает, что мы можем вычислить a, b и d точно. На практике же, вычисление a может выйти за пределы машинно-представимых чисел с плавающей точкой и, таким образом, мы получим только приближение a¢ к оптимальному наклону a. Затем мы должны выбрать b так, чтобы минимизировать максимум |a¢x
+ b - | вместо |ax
+ b - |. И опять мы сможем только вычислить приближение b¢ к оптимальному b. Ошибка аппроксимации d тогда будет максимумом |a¢x
+ b¢ - |; и снова мы способны вычислить только верхнюю границу d¢ для нее.

После этого формулы для расчета z
i
должны быть изменены для использования a¢, b¢ и d¢ вместо a, b и d. В действительности, вычисление z
0
, z
1
,…, z
n
по этим формулам будет ухудшено ошибкой округления; поэтому мы должны будем определить верхние границы d¢0
, d¢1
,…, d¢n
для этих ошибок и добавить их к ошибке аппроксимации d¢, всегда округляемой вверх, для получения ошибочного члена z
k
.

Аффинная оценка функции
e

x

Найдем линейную аппроксимацию ¦a
(e1
,…,en
) = z
0
+ z
1
e1
+ …+ z
n
en
для функции ¦(x
) = ex
. Пусть a, b и d имеют тот же смысл, что и в предыдущей части, [] = [a
, b
]. Тогда из (5) следует, что a = . Для вычисления b нам нужно знать максимум и минимум функции d
(x
) = ex
- ax
на [a
, b
]. Производная d
(x
) равна d
¢
(x
) = ex
– a. Следовательно, d
¢
(x
) = 0 при x
* = ln
(a)Î[a
, b
]. Но d
¢¢
(x*
) = ex
> 0, то есть x*
есть точка, в которой достигается минимум функции d
(x
) на [a
, b
]. Очевидно, что максимум достигается в точках a
или b
(причем из определения a имеем d
(a
) = d
(b
)). Таким образом, b = (d
(a
) + d
(ln(a)))/2 и d = (d
(a
) – d
(ln(a)))/2.

Как и в случае квадратного корня, мы сможем лишь получить верхние и нижние границы для величин a, b и d (здесь можно воспользоваться механизмом интервальной арифметики). После этого остается только немного перестроить формулы, по которым вычисляются z
i
, чтобы получить ¦a
.

Представление аффинных форм на ЭВМ

Описанный выше математический аппарат был применен автором в разработке конкретных библиотек на языке программирования C++. В них описаны такие классы, как интервал (interval), аффинная форма (affine_form) и смешанная аффинно-интервальная форма (mixed_form). Класс affine_form представляет аффинную форму , зависящую от m
символов шума. Он содержит такие поля как количество символов шума в форме, центр формы (т.е. x
0
) и массив из m
членов,
состоящих из частичного отклонения x
i
и соответствующего индекса i
– целого числа, которое уникально определяет символ ei
. Обозначения символа шума требуется записывать потому, что аффинная форма достаточно разрежена: поскольку данная программа может генерировать миллиард символов шума, каждая аффинная форма будет обычно зависеть от малого их подмножества. Поэтому ясно, что имеет смысл для каждой аффинной формы записывать только ненулевые члены x
i
ei
.

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

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

Общие подформулы

Часто вычисляемые выражения могут содержать некоторые подформулы более одного раза. Оценка таких выражений порождает многократное вычисление одной и той же подформулы, что является нерациональным. Кроме того, в AA
это может повлиять еще и на точность результата. Причина состоит в том, что вычисляемая в очередной раз оценка повторяющейся подформулы представляет ошибки ее линеаризации множеством новых символов шума. Таким образом, эти подформулы подразумеваются лишь частично зависимыми друг от друга, в то время как они тождественны. Это несоответствие предопределяет ошибки, которые могли бы быть отброшены на последующих шагах. Поэтому, когда кодируются выражения вида для аффинного оценивания, вдвойне важно определить общие подвыражения, подобно x
2
, и вычислить каждое из них только однажды.

Практическая реализация алгоритмов

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

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

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

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

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

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

Проблемы 1 – 4 так или иначе решаются путем применения интервальной или аффинной арифметики, методов, разработанных для автоматического анализа функций и с помощью методов структурирования геометрических объектов (подробности смотрите в следующих параграфах). Проблема 5 даже при поверхностном рассмотрении оказывается очень сложной и требующей дополнительного исследования, выходящего за рамки данной работы.

Автором разработаны соответствующие алгоритмы на языке программирования C++ для решения проблем 1 – 4, причем главная цель, которая здесь преследовалась, состояла не в реализации отдельных методов, а в построении общей вычислительной среды или системы, основанной на самых современных принципах программирования, в которой разработка конкретных методов оказалась бы тривиальной задачей. В качестве примера, демонстрирующего возможности этой системы, с помощью нее разработан алгоритм, реализующий метод ветвей и границ с градиентной проверкой, и протестирован на нескольких тестовых задачах.

Автоматический анализ функций

Пусть целевая функция подается на вход оптимизатора в виде символьной строки и списка имен переменных, от которых она зависит (для простоты здесь полагается, что функция задана в виде суперпозиции известных элементарных функций, а не алгоритмически). Если полученные данные корректны, то, производя лексический и синтаксический разбор входной символьной строки эту функцию всегда можно представить в виде дерева, с узлами, содержащими данные трех типов: константы (t_constant_data), переменные (t_variable_data) и операторы (t_operator_data). Как известно, между обычными и бинарными деревьями существует взаимно однозначное соответствие. Поэтому обычные деревья удобнее обрабатывать и хранить именно в виде право-прошитых бинарных деревьев [15], которые во всех дальнейших рассуждениях, для краткости, называются просто деревьями. На C++ такое дерево может быть описано в виде класса t_function_node, определив операции над объектами которого (по сути операции над деревьями) можно строить сколь угодно сложные суперпозиции функций. Например, если T
1
и T
2
деревья (или объекты класса t_function_node), то T
1
+ T
2
есть дерево с корнем, содержащем данные типа “binary_plus_operator” и имеющее в качестве левого поддерева дерево T
1
, а в качестве правого подде рева дерево T
2
(см. рис. 2).

cos

+

=

T 1 + T 2

+

/

y

x

ln

5

y

T 2

y

/

5

ln

y

x

T 1

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

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

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

Эффективное распределение памяти

При тестировании метода ветвей и границ на некоторых задачах обнаружилась интересная особенность. На первых этапах выполнения алгоритма скорость вычислений, произведенных с аффинными формами и с интервалами, в среднем равнялась нескольких тысячам в секунду. После того, как количество параллелепипедов в допустимом множестве достигало определенного количества, скорость вычислений неожиданно падала до нескольких сотен или даже десятков операций в секунду, даже при наличии свободной оперативной памяти. Этот эффект можно объяснить только тем, что алгоритму очень большое число раз приходится хаотично распределять и освобождать память под объекты небольшого размера, такие как параллелепипеды или вектора из аффинных форм. В результате оперативная память оказывается сильно сегментированной, и стандартный распределитель памяти становится недопустимо медлителен. Для решения этой проблемы был разработан собственный распределить памяти, базирующийся на стандартном. Общий принцип работы этого распределителя следующий. Память для объектов фиксированного размера size
распределяется заранее и представляется в виде списка сегментов некоторого, достаточно большого размера, которые в свою очередь разбиты на блоки размера size
, представленные в виде двух списков – свободного и занятого. Теперь, если программа затребовала кусок памяти размером size
, то распределитель просто переносит блок из свободного списка в занятый и возвращает на него указатель. Это происходит, конечно, только в том случае, если свободный список не пуст, иначе распределятся новый сегмент размером, равным сумме размеров уже распределенных сегментов (т.е. итоговый размер удваивается), разбивается на блоки размера size
которые, в свою очередь, добавляются в свободный список, после чего процедура повторяется. Механизм освобождения памяти происходит аналогично, нужно только определить момент, когда список сегментов не будет содержать занятых блоков и тогда освободить весь этот список.

На практике данный буферизированный распределитель оказался очень эффективным и легко применимым. Например, в одной из задач результат выполнения алгоритма содержал около полумиллиона параллелепипедов и функция разбиения этого множества на односвязные подмножества с помощью BSP-дерева (см. ниже), основанная на стандартном распределителе памяти затрачивала около 20 минут. При использовании же буферизированного распределителя потребовалось всего 6 секунд.

Структурирование результатов

На этапе завершения метода ветвей и границ результирующий список параллелепипедов может очень большим, поэтому получение такой информации как количество и содержимое связных областей, расстояние между этими областями, плотность распределения параллелепипедов в пространстве, границы результата и т.д. будет сопряжено со значительными трудностями. Поэтому нужен некоторый способ структурирования результирующего списка. В случае выбора двоичной схемы разбиения список удобно представлять в виде BSP-дерева, а в случае разбиения на k
частей, в виде k
-дерева. Такое дерево содержит узлы с информацией о параметрах разбиения и листья, являющиеся параллелепипедами из списка, и строится по мере выполнения алгоритма. Здесь возникает проблема построения сложных структур данных, когда один объект может одновременно принадлежать нескольким структурам данных (как, например, параллелепипеды принадлежат одновременно списку и дереву). Эта проблема может быть легко и элегантно решена путем применения механизма шаблонов языка C++.

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

Все описанные выше методы, а именно автоматический анализ функций, буферизированное распределение памяти и структурирование результатов были разработаны автором в виде библиотек на языке C++ и успешно протестированы на широком классе задач.

Численные результаты

Для сравнения возможностей интервальной и аффинной арифметики было проведено множество тестовых оценок примерных функций (некоторые из них показаны на графиках в приложении). Эти тесты показали, что способность аффинной арифметики учитывать взаимосвязи операндов обеспечивает в большинстве случаев значительно более точные оценки, чем в интервальной арифметике.

Ниже сравнивается выполнение IA
, AA
и AAIA
при решении некоторых хорошо известных задач глобальной оптимизации тремя видами интервальных BAB
методов: чистые модели, ускорение с локальным поиском; с ускорениями первого порядка (интервальная градиентная проверка); а также с Back-Boxing
’ом. Данные для таблиц, а также рисунки разбиения областей взяты из [7].

Поиски останавливались, когда в главном списке L не оставалось прямоугольников большего диаметра, чем определенная точность. При решении задач использовалась возрастающая точность (10-3
, 10-6
и 10-9
) для определения влияния точности на выполнение. В таблицах ниже приведено CPU время в секундах, необходимое для решения каждой задачи каждым из вариантов, а также число прямоугольников, остающихся в L после завершения. Время ¥ означает, что программа не была завершена в течении 15 минут CPU времени.

Также даются рисунки разложений области при использовании градиентного теста и точности 0.0625% от диаметра области. На них прямоугольники, показанные белыми, были отброшены, а серые прямоугольники остались после завершения.

Задачи были запущены на машине 166Hz Pentium под Linux, используя компилятор GNU g++.

1.
Функция Буфа
(таблица 1, рис.3)

¦(x
, y
) = (x
+ 2y
- 7)2
+ (2x
+ y
- 5)2

Глобальный минимум ¦ в квадрате W = [-10, 10]´[-10,10] равен ¦* = 0 = ¦(1, 3). Таблица 1 содержит результаты решения этой задачи всеми вариантами.

Эта задача была выбрана потому, что она выделяет главную слабость AA
: как обсуждалось ранее, AA
может выполняться плохо на функциях, которые являются суммами квадратов или когда зависимость между операндами незначительна. Это привело к большему времени для AA
, но только на множитель k
[2,3].

2.
Функция
Exp2
(таблица 2, рис. 4)

¦(x,y) = exy
(4x
2
+ 2y
2
+ 4xy
+ 2y
+ 1).

Глобальный минимум ¦ в квадрате W = [-10, 10]´[-10,10] есть ¦* = 0 = ¦(0.5, -1). Квадратичная часть этой функции содержит зависимости, которые позволяют AA
делать намного лучшие оценки. IA-
оценки возле глобального минимума оказались мало точны, что явилось следствием очень большого числа прямоугольников в списке для чистых BAB
методов. С другой стороны, AA-
оценки для экспоненциальной части были хуже, чем IA-
оценки; как следствие метод не смог отбросить много прямоугольников, удаленных от глобального минимума.

Таблица 1. Результаты выполнения алгоритмов для функции Буфа

BAB -метод		Чистый		С град. проверкой		С Back-Boxing ’ом
Точность	Ариф-метика	Время	колич.	Время	колич.	время	колич.
10-3	IA AA AAIA	0.07 0.18 0.12	4 10 4	0.07 0.16 0.11	4 5 4	0.00 0.01 0.02	1 1 1
10-6	IA AA AAIA	0.12 0.34 0.19	4 10 4	0.12 0.27 0.21	4 5 4	0.00 0.01 0.00	1 1 1
10-9	IA AA AAIA	0.21 0.52 0.29	4 10 4	0.20 0.43 0.32	4 5 4	0.00 0.00 0.00	1 1 1

Рис.3. Разбиение области при минимизации функции Буфа в
IA (слева
), AA (в центре) и
AAIA (справа).

Таблица 2. Результаты выполнения алгоритмов для функции
Exp2

BAB -метод		Чистый		С град. проверкой		С Back-Boxing ’ом
Точность	Ариф-метика	Время	колич.	Время	колич.	время	колич.
10-3	IA AA AAIA	496.07 23.50 0.41	20610 14 14	0.29 39.34 0.48	28 5 5	0.30 38.69 0.59	30 4 5
10-6	IA AA AAIA	¥ 23.43 0.69	36727 14 14	0.44 39.44 0.75	28 5 5	0.46 38.91 0.83	28 4 5
10-9	IA AA AAIA	¥ 36.61 12.25	36635 4922 4157	0.59 40.51 1.12	32 8 8	0.61 39.51 1.18	28 7 7

Рис.4. Разбиение области при минимизации функции
Exp2 в
IA (слева
), AA (в центре) и
AAIA (справа).
Заключение

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

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

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

Литература

1. Iwaarden R. Van, An Improved Unconstrained Gl
obal
Optimization Algorithm,
PhD
the­sis, Department of Mathematics, University of Colorado at Denver, 1996. URL h
ttp:
//ww
.cs.hope.edu/
~rvani
waa/phd.
ps.gz.

2. Moore R. E., Interval analysis.
New York, Prentice-Hall, 1966.

3. Glassner A., Space subdivision for fast ray tracing.
// IEEE Computer Graphics & Application, 1994.

4. IEEE Standard for Binary Floating-Point Arithmetic.
ANSI/IEEE Std
754-1985, IEEE, New York, 1985.

5. Comba J.L.D., Stolfi J.,
Affine arithmetic and its appl
ications to computer graphics,
in Proceedings of VI SIBGRAPI
(Brazilian Symposium on Computer Graphics and Image Processing), 1993, pp. 9-18. URL ht
tp://d
cc.unic
am
p.b
r/~stolf
i/.

6. Figueiredo L.H. DE, Surface intersection using affine arithmetic,
in Pro­ceedings of Graphics Interface '
96, May 1996, pp. 168-175. URL ftp:
//c
sg.u
waterloo
.ca/pub
/lhf
/gi
96.
ps.gz.

7. Figueiredo L.H. DE, Stolfi J., Adaptive enumeration of implicit surfaces with affine arithmetic,
Computer Graphics Forum, 15 (1996), pp. 287-296.

8. Figueiredo L.H. DE, Stolfi J., Iwaarden R. Van, Branch-and-bounds methods for unconstrained global optimization with affine arithmetic.

9. Farouki R., Rajan V.,
Algorithms for polynomials in Bernstein form,
Computer Aided Geometric Design, 5 (1988), pp. 1-26.

10. Gill P.E., Murray W., Wright H., Practical Optimization,
Academic Press, San Diego, CA,
1981.

11. Hansen E.,
A generalized interval arithmetic,
in Interval Mathematics, K.
Nickel, ed.,
no. 29 in Lecture Notes in Computer Science, Springer Verlag, 1975, pp. 7-8

12. Назаренко Т.И., Марченко Л.В., Введение в интервальные методы вычислительной математики
, 1982

13. Калмыков С.А., Шокин Ю.И., Юлдашев З.Х., Методы интервального анализа
, 1986

14. Шокин Ю.И., Интервальный анализ
, 1981

15. Кнут Д., Искусство программирования для ЭВМ, Т.1,
1976

16. Коллатц Л., Крабс В., Теория приближений. Чебышевские приближения и их приложения
, Москва “Наука” 1978

Приложение
Примеры оценивания диапазонов функций



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

Рис. 1. , ¦(¦(x
)) оценивается на отрезке [-1,1] c шагом 0.0125 в интервальной арифметике.

Рис. 2. , ¦(¦(x
)) оценивается на отрезке [-1,1] c шагом 0.0125 в аффинной арифметике.

Рис. 3.
, ¦(¦(x
)) оценивается на отрезке [-1,1.2] c шагом 0.0126 в интервальной арифметике.

Рис. 4. , ¦(¦(x
)) оценивается на отрезке [-1,1.2] c шагом 0.0126 в аффинной арифметике.