CC BY
69
УДК 681.5.011 +519.853.4 ББК 32.965 + 22.18
АТОМНАЯ ОПТИМИЗАЦИЯ, ЧАСТЬ 1: ТРАНСФОРМАЦИЯ ПРОСТРАНСТВА ПОИСКА И ОДНОМЕРНЫЕ ЗАДАЧИ
Поздяев В. В.1
(Арзамасский политехнический институт (филиал) Нижегородского государственного технического университета им. Р. Е. Алексеева, Арзамас)
Рассмотрены задачи оптимизации с полиномиальными целевой функцией и ограничениями в виде неравенств. Представлена трансформация основанного на теории моментов метода их решения, позволяющая конструировать эквивалентные алгоритмы решения в расширенном исходном пространстве поиска вместо пространства моментов. Детально рассмотрен случай одномерных задач оптимизации.
Ключевые слова: нелинейное программирование, матричные неравенства, полиномиальные неравенства, теория моментов.
Введение
Рассмотрим задачу нахождения глобальных экстремумов полиномиальной целевой функции на множестве, заданном полиномиальными неравенствами (ПН):
/ * = тіп / (х),
X
(1) 9г(х) ^ 0,
х є Мга, і = 1,..., т, где /(х) и 9і(х) — (не обязательно выпуклые) полиномы. Далее систему таких неравенств 9і(х) и саму данную задачу мы будем называть системой и задачей ПН соответственно; множество х,
1 Владимир Васильевич Поздяев, кандидат физико-математических наук, старший преподаватель (vpozdyayev@gmail.com).
удовлетворящих всем дг(х) ^ 0, — областью поиска или допустимой областью.
В теории управления большую роль играет аналог задачи ПН с участием полиномиальных матричных неравенств (ПМН, а также система и задача ПМН соответственно):
/ * = шш / (х),
(2) Сг(х) > 0,
х е Мга, г = 1,..., т, где Сг(х) — матрицы, элементы которых являются полиномами от х, а знак неравенства понимается как требование положительной полуопределенности. Очевидно, что можно рассматривать ПН как частный случай ПМН, а ПМН сводить к системе ПН с помощью, например, критерия Сильвестра. Следовательно, данные представления являются взаимозаменяемыми. Тем не менее, если при формализации задачи естественным образом возникают ПМН, то работа с ними напрямую (без сведения к ПН) может существенно упростить вычисления.
Отметим, что в литературе зачастую рассматривается не указанная общая форма ПМН, а ориентированные на те или иные классы задач теории управления конкретные их разновидности, такие как билинейные матричные неравенства [10] и параметризованные линейные матричные неравенства [1].
Работа с системами ПН и ПМН представляет собой трудную вычислительную задачу, известные алгоритмы решения которой неизбежно обладают такими недостатками, как ограниченность сферы применения, малая эффективность, эвристичность и т.п. В отдельных частных случаях (таких как задача о статической обратной связи по состоянию, динамической обратной связи по выходу без ограничения на порядок регулятора) были найдены подходы, позволяющие сделать задачу выпуклой и найти эквивалентные формулировки в терминах систем линейных матричных неравенств (ЛМН); см., например, [3, 11]. Для других, не менее фундаментальных, задач, таких как проектирование ПИД-регулятора, одновременная стабилизация, стабилизация статической обратной связью по выходу, эквивалентные системы ЛМН
неизвестны.
Задачи ПН и ПМН могут решаться локальными или глобальными методами.
• Локальные методы. Чувствительны к выбору начального приближения и не гарантируют сходимость в ряде важных случаев. В литературе описаны — и реализованы в виде программного обеспечения — различные варианты локальных методов для разных типов задач. В частности, первой общедоступной реализацией решателя билинейных матричных неравенств является пакет PENBMI [10, 13, 12].
• Глобальные методы. Зачастую основаны на методе ветвей и границ, но есть и иные варианты; в данной статье рассматривается подход на базе проблемы моментов. Требуют значительных вычислительных ресурсов. Важную роль здесь играет представление отдельных этапов решения задачи в виде ЛМН. Его стратегия может сильно влиять на эффективность метода, но, к сожалению, в общем случае здесь невозможно избежать той или иной разновидности комбинаторного взрыва. Как следствие, применимость данных методов ограничена задачами относительно небольшого размера.
В серии работ [7, 8, 14] был предложен и реализован в виде пакета для MATLAB [9] метод глобальной оптимизации, фундамент которого образуют теория разложения полиномов в сумму квадратов и двойственная ей теория моментов. С их помощью для задач ПН и ПМН строится иерархия ЛМН-релаксаций в пространстве моментов переменных. Каждый следующий элемент иерархии имеет более высокую размерность, но при этом гарантирует более точную оценку величины экстремума. Более того, теоретическая гарантия асимптотической сходимости к глобальному экстремуму на практике зачастую превращается в сходимость за конечное число итераций2. Для таких случаев предусмотрены
2 Итераций верхнего уровня, каждая из которых заключается в построении и решении ЛМН-релаксации определенного порядка.
процедуры, позволяющие извлечь положение глобальных экстремумов из решений ЛМН-релаксаций.
Отметим следующие достоинства данного метода:
• Позволяет решать задачи со скалярными или матричными полиномиальными неравествами.
• Гарантирует, что найденный экстремум действительно является глобальным.
• При наличии нескольких глобальных экстремумов находит их всех одновременно без привлечения дополнительных стратегий перебора, ветвления и т.п.
• Основная масса вычислений связана с решением систем ЛМН определенного вида, для чего могут быть использованы эффективные реализации известных алгоритмов.
Вместе с тем он имеет и существенные недостатки, в том числе:
• Он также не избежал комбинаторного взрыва: например, в [8] приводится пример задачи о стабилизации статической обратной связью по выходу линейной системы с 9 координатами состояния, 1 входом и 5 выходами, для которой попытка решения на основе неравенства Ляпунова привела бы к ЛМН с 33649 неизвестными уже в простейшей (первой) релаксации.
• Максимальный порядок моментов, задействованных в ЛМН-релаксации, растет вместе с ее номером, равно как и со степенью полиномов в исходной системе ПМН. При этом ухудшается обусловленность матриц, составляющих ЛМН-релаксации, что отрицательно влияет на точность решения. Эта проблема дополнительно усугубляется потенциально плохой обусловленностью самих ПМН; см., например, [6].
Автором была разработана трансформация данного метода, возвращающая нас из пространства моментов неизвестных (используемого в ЛМН-релаксациях) обратно в исходное пространство поиска (используемое в ПМН) с некоторыми дополнениями. Такой переход позволяет устранить оба указанных недостатка, сохранив при этом способность находить как локальные, так и глобальные экстремумы в задачах с невыпуклыми и даже несвязными областями поиска.
Отметим, что многие подвиды задачи ПМН, а следовательно, и она сама, являются ЛГ’-трудными [2]. Это не позволяет нам сконструировать алгоритм полиномиальной сложности, гарантирующий нахождение глобального минимума в общем случае. С другой стороны, если рассматривать не сформированные специальным образом, а возникающие на практике задачи, можно видеть, что даже в случае их невыпуклости рельеф целевой функции и вид области поиска могут оказаться относительно простыми (см. иллюстрацию к разделу в [8]). Вышеуказанный метод лишен возможности воспользоваться этим обстоятельством, но предлагаемый в данной работе алгоритм имеет достаточную гибкость, чтобы, варьируя количество задействованных атомов, балансировать между объемом вычислений и результатами, характерными для локального поиска и для полноценного эквивалента исходных ЛМН-релаксаций (см. разделы 4 и 5).
В настоящей работе детально рассматривается случай скалярных неравенств и одномерного пространства поиска: задача ПН с п = 1. Несмотря на очевидную ограниченность, он позволяет полностью сформировать структуру нового алгоритма оптимизации, а также представить базовые формулы, получающие свое развитие в задачах с п > 1, в том числе и в ПМН общего вида.
В разделе 1 кратко описан исходный метод построения ЛМН-релаксаций для задач ПН. Раздел 2 содержит результаты, связанные с трансформацией пространства поиска данного метода. В разделах 3 и 4 представлено применение этих результатов к случаю п = 1, дающее вычислительную схему, значительно бо-
лее гибкую по сравнению с оригиналом и при этом имеющую аналогичные характеристики сходимости. Пример, сравнивающий базовый метод с его трансформацией, приведен в разделе 5.
В дальнейших работах схема, описанная в разделах 2-4, будет распространена на задачи ПН с n > l, а также на задачи ПМН.
1. Известные результаты
Приведем основные положения метода решения задач ПН, опубликованного в [14] и расширенного в [7, 8].
Введем следующие обозначения. Пусть br (x) — вектор, состоящий из одночленов, образующих базис пространства многочленов порядка не выше r:
br (x) = [ l x1 ... xn x2 x1x2 ... x" ... x" ]T ,
а sn(r) = Cn+r = ^ — его размерность. Каждому одночлену из br (x) поставим в соответствие вектор а є N", ^ i ai ^ r (далее будем записывать как a ^ r), показателей степеней x1, x2, ..., x^ обозначим xa = x"1 x"2 .. .x"n. Для произвольного многочлена p(x) степени не выше r будем индексировать вектор p его координат в базисе br (x) двумя взаимозаменяемыми способами: по номеру элемента и по вектору показателей степеней. Таким образом, p = bi]1<i<sn(r) = [PajaGNg, a^r, в том числе P1 = P[0,0,...,0], P2 = P[1,0,...,0] и т. д., при этом
p(x) = (p, br(x)) = pTbr(x) = ^Paxa.
a^r
Аналогичным образом будем индексировать строки и столбцы матриц там, где это применимо.
Рассмотрим вектор y = [ya]aeNg, a^2k и полином 9(x) степени не выше 2d, d є N0, d ^ k. Построим матрицу моментов Mk(y) є RSn(k)xSn(k) и локализующую матрицу Mk-d(9,y) є Msn(k-d)xsn(k-d), строки и столбцы которых проиндексированы
соответственно базисам Ьд(ж) и (ж), а элементы заданы сле-
дующими формулами:
(М (У))ав = У«+в ,
(3) (Мк_^(д,У))ав = 5^ ^У«+в+7,
где а, в, 7 ^ МП, а д7 — элементы вектора коэффициентов д(ж) в базисе Ь2^(ж).
Нас интересуют ситуации, когда для у существует положительная Борелевская представляющая мера ^ такая, что
уа = J жа ё^, а є МП, а ^ г.
Величины уа при этом являются моментами ж, а сами векторы у принадлежат (усеченному) пространству моментов ж. Матрицы
Мд, (у) и Мк_й (д, у) тогда равны
Мй (у) = Ь (ж)Ьк (ж)т ф,
(4) 7,
Мк_^(д,у) = д(ж)Ьд_ (ж)^_гі(ж)т ф,
и мы можем сформировать два условия, играющие центральную роль в дальнейших построениях.
• Необходимым условием существования меры с моментами у является Мд(у) ^ 0. В самом деле, для произвольного многочлена р(ж) степени не более к (так что р — произвольный вектор из М5п(д)) имеем:
/к*)2 Ф =
= У (рА(ж))2 ф = J (рТЬй(ж)) (^(ж)тр) ф =
= У рт (Ьд(ж)Ьй(ж)т) рф = Р^У Ьд(ж)Ьй(ж)т ф^ р =
= ртмй (у)р.
Поскольку левая часть неотрицательна, правая должна представлять собой положительно полуопределенную квадратичную форму.
• Аналогично можно продемонстрировать, что для того чтобы носитель ^ был подмножеством {ж € Rn|g(x) ^ 0}, необходимо выполнение неравенства Mfc_d($, у) ^ 0:
y"g(x)p(x)2 = Jg(x) (p,bfc-d(x))2 =
= J g(x) (pTbfc-d(x)) (bfc-d(x)Tp) d^ =
= Р^У g(x)bfc-d(x)bfc-d(x)T d^ p =
= pTMfc_d(g,y)p ^ 0.
Здесь p(x) — произвольный многочлен степени не более k — d (так что p — произвольный вектор из RSn(k_d)).
Вернемся к задаче (1). Пусть dj = [| deggj(x)"|, а k удовлетворяет ограничениям 2k ^ deg f (ж), k ^ di. Пусть fa — коэффициенты f (ж) в базисе b2k(ж), так что
f f (x)d^ = f V fax“ d^ = V faya.
J J a^2fc a^2fc
ЛМН-релаксацией (1) будем называть систему
f * = min V f«y«, y
a^2fc
(5) Mfc(y) ^ 0,
Mk-di(gi,y) ^ 0, i = l,...,m, y[0,0,...,0] = 1.
В [14] (теорема 4.2) показано, что, с учетом некоторых непринципиальных ограничений, при k ^ то величина экстремума ЛМН-релаксации стремится к величине экстремума исходной задачи
ПН. Более того, как правило, уже при конечных (и относительно небольших) значениях k данные величины становятся равны, а вектор моментов решения задачи ПН является решением соответствующей ЛМН-релаксации.
Согласно [З] (теорема 1.6), достаточным условием достижения вышеуказанного значения k является
r = rank Mk (y*) = rank Mk-d(y*),
где y* — решение ЛМН-релаксации, а d = maxi di. Если оно выполняется, то у* представляет собой вектор моментов r-атомной меры3, атомы которой x*j, j = l,..., r, соответствуют глобальным минимумам (1). Данные атомы могут быть найдены из системы полиномиальных уравнений, следующей из (4):
Г
Mk(y*) = ^ Ajbk(x*j)bk(x*j)T, j=i
где Aj — неизвестные веса, удовлетворящие условиям Aj ^ О, j Aj = l. В [7] представлен алгоритм эффективного решения таких систем с помощью базовых операций линейной алгебры.
В [8] показан способ построения аналогичных ЛМН-релаксаций для задач ПМН без предварительного перевода в форму ПН.
2. Трансформация пространства поиска
Изложенный в предыдущем разделе метод решения задач ПН включает три этапа (первые два шага могут повторяться, пока мы не найдем подходящее значение k).
1) Прямая трансформация пространства поиска из R" в RSn(2k): построение ЛМН-релаксации (З) с достаточно большим k.
3 N-атомная мера — мера, носитель которой является множеством из N точек («атомов»).
2) Решение вспомогательной задачи (полученной системы ЛМН).
3) Псевдообратная трансформация пространства поиска из ^п(2к) в ^гх(п+1): восстановление экстремумов исходной задачи из полученного решения.
Таким образом, в решении задачи участвуют три различных пространства:
• исходное пространство поиска (Мга);
• пространство моментов порядка до 2к (М5п(2к));
• пространство атомов, элементами которого являются кортежи из г пар «вектор исходного пространства; его вес»
(Мгх(га+1)).
Выбор пространства моментов в качестве целевого пространства для этапа трансформации позволяет линеаризовать задачу оптимизации, но приводит к возникновению существенных затруднений нового характера:
• Размерность пространства моментов 8га(2к) подвержена комбинаторному взрыву с ростом п и к.
• При увеличении к нередко возникает тенденция к появлению существенной разницы в порядках величин уа, связанная с возведением элементов х в высокие степени. Это отрицательно сказывается на точности вычислений как на этапе решения ЛМН-релаксации, так и при извлечении атомов из найденного вектора моментов.
Чтобы справиться с ними, в настоящей работе мы предлагаем исключить из процедуры решения выход в данное промежуточное пространство. Для этого необходимо сформулировать (и решить) напрямую в пространстве атомов вспомогательную задачу оптимизации, в определенном смысле эквивалентную ЛМН-релаксации (5). При этом мы отказываемся от ЛМН-48
представления исходной задачи, заменяя его на структуру, объединяющую непосредственно элементы (1) в рамках одной целевой функции, минимум которой находится с помощью модифицированного метода Ньютона.
2.1. АЛГОРИТМЫ ОПТИМИЗАЦИИ
Прежде чем формулировать задачу оптимизации, эквивалентную (5), определим, что подразумевается под «эквивалентностью» в данном контексте. Запишем систему ЛМН (5) в общем виде:
/*& т
/ = шт с у,
у
(6) Р(у) = ^ Р«у« ^ 0
у[0,0,...,0] = 1
где вектор с состоит из /а, а Р (у) =
diag(Mfc(у),Мй-гіі(дьу),... ,М^т(дт,у)). Для ее решения существует множество алгоритмов, см. [4]. Далее мы будем использовать один из простейших вариантов метода внутренней точки в прямой форме с ньютоновским направлением поиска. Приведем его краткую схему.
Выберем4 вектор у(0), удовлетворяющий ограничениям (6), и монотонно невозрастающую сходящуюся к 0 вещественную последовательность {^(г)}, і = 0,1, 2,.... Сформируем семейство комбинаций целевой и барьерной функций вида
/(і) (у) = сту - ^(г) ^ det Р(у).
Для і = 0,1, 2,... положим у(г+1) равным результату решения подзадачи поиска минимума /(і) (у) на интервале прямой (с учетом ограничения Р(у) ^ 0), проходящей через у(і) и имеющей направление, предписываемое методом Ньютона. Конструируемая таким образом последовательность {у(і)} сходится к искомому решению задачи (6).
4 Мы не останавливаемся на вопросах выбора начального приближения и стратегии формирования {м(г)}> поскольку они подробно изучены в литературе.
Подзадача о линейном поиске представляет собой задачу одномерной оптимизации функции вида /(г)(у(г) + £ Ду(у(г))) относительно £, где Ду = Ду(у) — направление спуска. Отметим, что, поскольку У[о,о,...,о] является константой, при формировании Ду необходимо использовать градиент и гессиан /(г)(у) в подпространстве М5п(2к), соответствующем всем уа для а = 0:
Ду = [ о | Ду ],
Ду = - (V2/«(у))- V/«(у),
где V/(г) и V2/(г) — градиент и гессиан / (г)(у) в указанном под-пространстве(: )
(V/(%)) = Са - ^ ^ (Р(у)-1Р«) ,
(7) а
(V2/«(у)) = ^ ^ (Р (у)-1ЫР(у)-1Рв) .
\ / ар
В силу неравенства Р(у + £ Ду) = Р(у) + £ Р(Ду) ^ 0, интервал поиска ограничивается сверху и снизу ближайшими к 0 положительным и отрицательным обобщенными собственными значениями матриц Р(у) и —Р(Ду) (если таковые имеются).
В качестве результата решения данной подзадачи может использоваться вектор у; = у + £*Ду, где £* — найденное положение экстремума (исчерпывающий поиск), или же у; = у + а£*Ду для некоторого а € (0; 1) (демпфированный поиск). Последний вариант в некоторых случаях позволяет улучшить сходимость и устойчивость алгоритма.
Пусть теперь у является гладкой функцией вектора х той же размерности5: у = у(х); матрица Якоби Т с элементами
Т = Йу*
= ёх,-
5 Впоследствии в роли х будет выступать вектор из пространства атомов.
невырождена, и (конечные) минимумы /(г) (у) и / (г)(у(х)) существуют. Для малых значений £ имеем:
(8)
уг(х+£ Дх) = уг(х)+^ Дхз£+о(£) = уг(х)+(Т Дх)г£+о(£).
3 3
Таким образом, движение от точки х в направлении Дх эквивалентно в малом движению от у(х) в направлении Ду = Т Дх.
Такое понятие эквивалентного направления естественным образом интегрируется в алгоритм оптимизации. А именно, построим алгоритм поиска экстремумов семейства функций / (г)(у(х)) в пространстве векторов х, повторяющий вышеприведенную схему, в котором
• начальным приближением является вектор х(0) такой, что
у(х(0)) = у(0);
• направление спуска выбирается как Дх = Дх(х) =
т-1Ду(у(х)).
Данный алгоритм мы будем называть эквивалентным исходному алгоритму в контексте трансформации у = у(х). Далее мы показываем, что, в дополнение к вычислению / (г)(у(х)) как функции от х, вектор у можно также исключить из функции Дх(х) = Т-1Ду(у(х)), а следовательно, и из нового алгоритма в целом.
Отдельно отметим связь между производными целевых функций по х и по у:
Й /(г) = / (г)
ёхг ^ ёхг ёу3 ,
з
или -дх /(г) = ТТ_ау /(г). Отсюда следует, что для найденного новым алгоритмом экстремума х* (в котором -дХ/(г) (у(х*)) = 0) вектор у* = у(х*) (для которого в силу невырожденности Т необходимо имеем -ду/(г)(у*) = 0) является решением задачи (6). Это верно для любого допустимого начального приближения — в том
числе и тогда, когда область поиска в исходной задаче (1) несвязна, а решение задачи и выбранное в качестве начального приближения множество атомов находятся в разных компонентах связности: см. раздел 5.
2.2. БАЗОВЫЕ РЕЗУЛЬТАТЫ
Прежде чем представлять формулу для нахождения Дх(х), приведем два вспомогательных результата.
Для целевой функции / (х), которая в окрестности некоторого хо имеет вид
/ (х) = (х — х0)ТА(х — х0) + ЬТ(х — х0) + с + о(||х — х0||2),
следующая лемма позволяет найти ньютоновское направление с учетом линейного ограничения в виде равенства.
Лемма 1. Пусть в задаче оптимизации
/ * = шш / (х),
/ (х) = (х — х0 )ТА(х — х0) + ЬТ(х — х0) + с,
ТТ
V х = V х0,
А € Мгахга; х,х0,Ь, V € Мга/ с € М, выполняются соотношения Ь, V € со1 А + АТ/ А + АТ ^ 0. Тогда вектор Дх, который необходимо прибавить к х0, чтобы попасть в точку минимума (одну из возможных, если det А + АТ = 0), равен
где д = V/(х0) = Ь, Н = V2/(х0) = А + АТ, Н — произвольная 6 обобщенная обратная к Н матрица.
6 Свободно выбирается из множества обобщенных обратных матриц при первом упоминании. Все дальнейшие вхождения этой матрицы имеют то же значение.
Доказательство. Воспользовавшись методом множителей Лагранжа, получаем систему уравнений
Н (ж — жо) + д + = 0,
т т
V тж = V тж0,
где х — искомая точка минимума, А — множитель Лагранжа. Поскольку д, V € со1 Н, первое уравнение разрешимо относительно
ж:
ж = ж0 — Н (д + Аv).
Подставив результат во второе уравнение, находим
А=
vTH д
V тН—V,
откуда
ж = жо + Н —д +
vTH д
--------—1
V тН — V
Учитывая, что Дх = х — х0, получаем утверждение леммы.
Следующая лемма показывает влияние трансформации у = у(х) на градиенты и гессианы функций.
Лемма 2. Пусть х € Мп, у € Мп ^ Мп, / : Мп ^ М, и функ-
ции у(х) и / (у) достаточное количество раз дифференцируемы. Тогда
Т
дх = Т ду,
нж = т Т Ну т + ну,
где
дх
дУ = /,
%
ёж
Ух/, Нх = V/ НУ )
Ну = уу /,
.7 J
V и V2 с нижними индексами — полные градиент и гессиан по соответствующим переменным:
=
ё
ёр,
ё2
ёр,
■9
у
и функции в правых частях вычисляются в точках х (для у(х), /(у(х)) и их производных) и у(х) (для /(у) и ее производных).
Доказательство. Данные равенства являются матричной формой соотношений
ёж,-
£
ё/ ёу, ёу, ёж,
ё2/ ёж,ёж,
_ё_
ёж.
ё/ ёук ёук ёж,
ё ё/
£
1 к
У, ,
^ ' ёж, ёук / ёж
£
ёуь ё/ ё2ук
+
ёук ёж, ёж,-
ё2/ ёу ёук
ёуг ёук ёж, ёж
£+£
ё/ ё2ук
ёук ёж, ёж,
(суммирование производится по всем допустимым значениям индексов).
Запишем решаемую методом внутренней точки подзадачу оптимизации в несколько более общей форме:
/ * = шш / (у),
У
(9)
Т Т (0)
^У = !,У( ;,
где / = /, а = [100 ... 0] (из (6) также следует, что
^У(0) = 1, но на данном этапе это не имеет значения). Тогда имеет место следующая теорема (мы продолжаем использовать обозначения из леммы 2).
и
Теорема 1. Пусть в точке у = у(х) выполняются соотношения ду,^у е со1 Ну; Ну ^ 0; det 3 = 0. Тогда направление поиска Ах в алгоритме, эквивалентном методу внутренней точки в первичной форме с ньютоновским направлением поиска в контексте трансформации у = у(х), применительно к (9) имеет следующий вид:
где = Н, — Н,; Н- — обобщенная обратная к матрица; V, = 3т^у; градиент и модифицированный гессиан Н, вычисляются в текущей точке х.
Доказательство. Воспользуемся леммой 1:
величины ду = 3 тдХ и Ну = 3 т(НХ — Ну)3 1 = 3 тНХ3 \
Здесь использован тот факт, что для произвольной обобщенной
(10)
Ах = Нх
Подставим в Ах = 3 1Ау выраженные из результатов леммы 2
а также ^у = 3 т^Х:
обратной к матрицы Н- матрица 33Н,- 3т удовлетворяет соотношению
(3-тН,3-1)(3зН-3 т)(3-тя*3 х1) =
= 3хтяХяхяХ3х1 = 3-тяХ3-1,
55
а следовательно, является обобщенной обратной для
З-ТЯхЗ-1 — и наоборот, для произвольной обобщенной обратной к З-ТЯхЗ-1 матрицы (З-ТЯх З-1)- найдется соответствующая обобщенная обратная к Ях матрица
Я- = З-1 (З-ТЯх З-1)-З-Т:
ЯхЯ-Ях = (ЗТЗ-ТЯх)(З-1(З-ТіїхЗ-1)-З-Т)(іїхЗ-1З) =
= зТ(з-ТЯхз-1)(з-ТЯхз-1)-(з-ТЯхз-1)з =
= ЗТ(З-ТЯхЗ-1)З = Ях.
В (10) есть только два элемента, которые содержат у или его производные в явном виде.
• Вектор ^х = ЗТ^У. Если у — вектор моментов, то, как мы увидим в следующем разделе, ^х почти столь же тривиален, как и .
• Матрица Яу = ^(Ухук). Чтобы исключить из данного выражения у, необходимо учесть фактический вид целевой функции:
/(,)(у) = сТу - ^(,) ^ё^ Р(у),
откуда
Яу = УУ / (,)(у) = -^(,)УУ logёet Р (у).
Таким образом, линейная компонента сТу не входит в Яу, а следовательно, и в Ях. Вычисление Яу для единственной имеющей значение компоненты, а именно, logёet Р(у), будем производить на основе следующей теоремы.
Теорема 2. Для функции / вида /(у) = logёet Р(у), где Р (у) линейно зависит от у, поправка в модифицированном гессиане Ях из теоремы 1 равна
ЯУ = £ ^ Т(ух^<,).
к i,j
где Ру = Р(у(ж))у и р_т = Р(у(ж))_т — соответствующие элементы матриц Р и Р_т (суммирование производится по всем допустимым значениям индексов).
Доказательство. Воспользовавшись известной формулой для матричной производной определителя (см., например, [15]), получаем:
¥- = ‘г (р-1 тг^ ) = £ Р_т 4^ •
ЙЩ; \ Гук/ 7^ у Гук
откуда
Г/ = у-л / у-л р_т
^ Йу; ёжр %3 Гук
= р-т ^у^ г2ук \
Т! ^ \к ГУк ГХР Гж«/
для всех допустимых р и д. Сравнивая последнее подвыражение в скобках с правой частью второй формулы из леммы 2 и учитывая, что УуРу = 0 (поскольку Р линейно зависит от у), видим, что
данное подвыражение — не что иное, как ^ ^, откуда
у^ Г/ ё2ук = р_т й2ру
^ ёук ёжр у ёжр
к
Доказываемая формула представляет собой матричную запись данного равенства.
3. Задачи одномерной оптимизации
Применим полученные в предыдущем разделе результаты к задаче (1) с п = 1 и ее ЛМН-релаксации (5), взяв в качестве трансформации пространства поиска переход из пространства атомов в пространство моментов. Элементы пространства
57
Й2 ук
ёжр
атомов будем обозначать как z, трансформацию — как y = y(z), а направление поиска в новом пространстве, соответственно, как Az.
Пусть z — вектор в пространстве атомов Rrx(n+1) = R2r, имеющий структуру
z = [ X1 X2 ... Xr P1 P2 ... Pr ] = [ XT I pT ]T ,
где xi Є R — атомы, pi є (0; І) — их веса. Вектор моментов у
имеет размерность S1(2k) = 2k + І и состоит из элементов вида
r
yj = У[,-1] = Y^ Pixj-1> j = І> 2,...>2k + І.
i=1
Вышеизложенные результаты неприменимы напрямую к векторам y и z, поскольку размерности последних не могут быть равны в силу различающейся четности. Чтобы исправить положение, необходимо или расширить теорему 1, распространив ее на случай отличающихся размерностей задействованных пространств, или модифицировать сами пространства. Наиболее адекватным вариантом нам представляется добавление еще одного момента — порядка 2k + І — в вектор у в сочетании с условием r = k + І. Новый момент не входит ни в (5), ни в Ay, и нужен исключительно для выполнения условий теоремы 1. Далее мы будем работать именно с таким расширенным вектором У є R2k+2.
Отметим, что (при n = І) для любой Mk (y) ^ 0 найдется соответствующая rank Mk (y)-атомная мера: см. алгоритм, изложенный в [7]. Поскольку rank Mk (y) ^ s1 (k) = k + І = r, каждая такая мера — а следовательно, каждая такая Mk (y) — может быть представлена некоторым вектором z (в котором, возможно, часть весов pi равна О).
Найдем матрицу Якоби J = J(z) трансформации y = y(z) =
у([хт|рт]т), а также вектор ^ = 3Т^ (здесь а = 2г — 1):
3
тТ
0 ... 0 1 ... 1
Р1 ... Рг х1 ... хг
1 аріж: а— 1 гу г) 'у -х ± • • • (_л г ^-/ г ™а X1 ... жа
1 0 ... 0 Н 0 . . 0 1
1 І'
Последний вектор содержит Г нулей и Г единиц. Значения Уу и ^ естественным образом отображаются на ограничения в соот-
ветствующих задачах оптимизации:
^ У
т
= 1
= 1,
г=1
Для применения теоремы 1 необходимо выполнение условий ду, уу Є сої Ну; Ну ^ 0; det 3 = 0. Поскольку добавленный нами момент У[2к+1] не входит в (5), а следовательно, и в /(і)(у), градиент и гессиан последней имеют структуру
$У
#У
я„
[ Н 0 '
0 0
где ду е М51(2к) и Ну е м«1(2к)Х51(2к) — градиент и гессиан / (г)(у), какими они были бы при использовании оригинального вектора моментов. Метод внутренней точки применительно к (5) обеспечивает в допустимой области соотношение Яу > 0, откуда получаем, что Яу ^ 0, а также что любой вектор из М81(2к+1) с нулевым последним элементом — каковыми являются ду и уу — принадлежит со1 Яу.
Формула для вычисления det 3(г) приведена в следующем утверждении. Сравнивая ее с det Мк(у(г)) из утверждения 3, видим, что det 3(г) =0 тогда и только тогда, когда det Мк(у(г)) = 0. Таким образом, внутри допустимой области det 3(г) = 0.
V
£
0
На ее границе данное неравенство может не выполняться, что способно привести к преждевременной сходимости алгоритма к неправильной точке, если последовательность приближений |^(г)} слишком быстро подходит к границе допустимой области. Последнее возможно при использовании последовательности излишне быстро сходящейся к 0. Таким образом, стра-
тегия формирования {^(г)} в новом алгоритме не менее важна, чем в исходном методе внутренней точки.
Введем дополнительные обозначения:
• для множества в и натурального числа п обозначим как
СП множество п-элементных подмножеств в; пусть также
пи пи
С[г] = С{1,2,...,г};
• для конечного числового множества в обозначим как в^ его г-й наименьший элемент: в1 — минимальный элемент в, в2 — следующий по величине и т. д.
Имеют место следующие утверждения.
Утверждение 1. Якобиан у(г) равен
( \
det 3(*) = (—1) ” Щ р* ) Л (ж*2 — ж*! )4
Уесм )
Чг=1
(для г = 1 последнее произведение считаем равным 1).
Доказательство. Для г = 1 определитель находится непо-
средственным вычислением. Предположим теперь, что г > 1, и доказываемая формула верна для г' = г — 1. Применим к 3 ряд элементарных преобразований. (В каждом пункте указанные операции применяются к строкам/столбцам матрицы, полученной в предыдущем пункте.)
• Для 3 = 2,..., г: вычтем из столбца 3 + г столбец г + 1.
• Для 3 = 2,..., г: вычтем из столбца 3 столбец 3 + г, умноженный на ———.
Xу — XI
• Для і = 2,..., г: вычтем из строки і строку і — 1, умноженную на ж1.
В результате приходим к матрице следующего вида (а = 2г — 2):
о 0 1 0 \
Р1 0
3' = Р1Ж1 31 0 32
у р1жа 0 /
где
31 =
I о
Р2(Ж2 — Ж1)
\ ар2Жа-1(Ж2 — Ж1)
32
/ Ж2 — Ж1 Ж2(Ж2 — Ж1)
о
\
Рг (жг — Ж1)
аргжа-1(жг — ж1) У
Жг — Ж1 \
Жг (жг — Ж1)
\ жа(ж2 — ж1) ... ж^(хг — ж1) у
Проделаем аналогичные манипуляции с целью упрощения не только первой, но и второй строки матрицы.
• Для 3 = 2, . . . , г: вычтем из столбца 3+ г столбец 1, умно-
женныи на
—Х1
Р1 ‘
• Повторим последние два пункта из предыдущего списка. Получаем финальную матрицу (а = 2г — 3):
0 0 1 0
Р1 0 —Ж1 0
0 0
0 31 0 32
0 0
где
О
Р2(Х2 - Жі)2
У ар2Ха 1 (Ж2 - Хі)2 .
7 ''
и о
( (Х2 - Жі)2
Х2(Х2 - Жі)2
У ж2(ж2 - Х1)2
О
Рт(Жт - Жі)2
артж2 і(жт - Хі)2 /
(жт - Хі)2 ^
Жт(жт - Жі)2
Ж2 (ж Х )2 і
і/ /
Разложив по первым двум строкам det 7'' и возникающий при этом минор, имеем (а = 2г — 3, х' = [ж2 ... хг р2 ... рг]Т):
det 7(г) = det 7'' = (- 1)Трі det ( 7і' | 72') = (- 1)Тріх
/О ... О 1 ... 1 \
Р2 ... Рт Ж2 . . . Жт
: Л (ж —жі)4 det
і=2
«а—і
—і /у>а
\ ар2ж2 ... арт ж а
(-1)тРі (II(Жі - жі)М 7(г')
= (-1)
і=2
г(г+1)
/
2
П рі
і=1
\
(жі2 - жі1)4
\іЄСИ У
поскольку
і=2
Л (жі2 - ЖІ1 )
\іЄС22,...,г} /
Прежде чем перейти к вычислению остальных элементов Ах, отметим, что матрица моментов Мк (у) может формально 62
рассматриваться как локализующая матрица М*_^0(до, у), где до (ж) = 1, 4 = 0. Кроме того, в силу структуры Р(у),
т
logdet Р (у) = logdet М* (у) + ^ log det М*_^ (д»,у),
»=1
так что
/(г)(у) = сту - ^(г) ^logdet М*(у) + ^ logdet М*._*(д*,у)^ ,
и величины и Н, из которых с помощью теоремы 1 находится Ах(х), равны значениям следующих функций в точке х:
(11)
т
/(г) = / - ^(г) ( logdet М* + ^ logdet М*_^(д»,у)
»=1
т
V2/М = -^(г) V2 logdet М* + ^ '72 logdet М*_^(д»,у)
2 -
»=1
где V2 = V2 — ^г^2у») -дут. Чтобы вычислить правые части, необходимо найти сами матрицы М*_^(д», у(х)), гессианы их элементов (в соответствии с теоремой 2), а также logdet М*_й.(д»,у) и V2 log det М*_й. (д»,у). Соответствующие формулы, не требующие вычисления вектора моментов у, представлены в следующих утверждениях.
Утверждение 2. Локализующая матрица Мк_^ =
М*-^(д,у(х)) равна
Г
М*_й = ^ р»д(ж»)Ь*_й(ж» )Ь*_й(ж»)т.
г=1
Доказательство. Следует из (4) и того, что ^ — г-атомная мера.
Утверждение 3. Определитель матрицы Мк_^ =
Мк_^(д,у(х)) имеет вид
det М*_й = ^ £>8,
где S = CjT] d, и
Ds = Щ P*g(x*) J I Л (xi2 - Xii)
Vies / \ieC2
Доказательство. Пусть BCT и BCT, где a С {1, 2,..., r}, — матрицы (r — d) x |a|, i-е столбцы которых равны bk—d(xCTi) и л/p^i)bfc-d(xCTi), соответственно. Несложно убедиться, что
— —T
^fc—d B{1,2,...,r}B {1,2,...,r}
и
BBCT = BCT diag (pCT1 g(xCT1 ),pCT2g(xCT2),...) BT Применяя формулу Бине-Коши, имеем:
(12) det Mk—d = det B{1,2,...,r}BT12,...,r} = ^ det B . det BT =
s€S
= ^ det B. | ^Pig(x) j det BSf.
s€S \i€s /
Поскольку (для n = 1) bk—d(®j) = [ 1 X ... xk—d ] , матрица B. является матрицей Вандермонда, и det B. = det BT = Пгес 2 (xi2 — xil ). Подставляя в (12), получаем доказываемую формулу.
Утверждение 4. С учетом обозначений предыдущего утверждения, градиент матрицы Mk—d = Mk—d(g,y(z)) может быть представлен как
Vz log det Mfc—d = gz—,
det Mfc—d
где
gz = Vz det Mfc—d = ^ (Vz log D.) D.,
s€S
а ненулевые элементы Vz log Ds имеют вид
-f log D. = , , , e ,.
ПО* • П ПГ ■ \ \ ^^
dxi g(xi) \ . ^ ж* - ж,-
\jes\{t}
d1 -—logD. = —, i e s. dpi Pi
Доказательство. Данные формулы получены непосредственным дифференцированием с использованием соотношений
f1
вида (log f )x = (применительно к f = log det Mk—d) и (exp f)X = (exp f)fX (применительно к f = log Ds).
Утверждение 5. С учетом обозначений предыдущего утверждения, гессиан матрицы Mk—d = Mk—d(g,y(z)) может быть представлен как
V2 log det Mk—d = Hz gzgT
det Mfc—d (det Mfc—d)2’
где
Hz = V2 det Mfc—d =
^ (V2 log D. + (Vz log D.)(Vz log D.)T) D.,
s€S
а ненулевые элементы V2 log D. имеют вид
d2 >ogd.=^2—i e
dxi dx- g(xi) \g(xiW \fcel\H (xi — Xk)2
{i,j} С s, i = j;
d2 2
log D. = 7-------------r,, {i,j }C s, i = j;
dxi dxj (xi — ж- )2
d2 1
— logD. = -pj, {i,j)C s, i = j.
Доказательство. Данные формулы получены непосредственным дифференцированием с использованием соотношений
f// f / f /
вида (log f )Xy = f------(применительно к f = log det Mk—d)
и (exP f )XXy = (exP f )(f"y + ff) (применительно к f = log D.).
2
Для матрицы моментов М& (у(я)) данные результаты существенно упрощаются, поскольку д0 (х) = 1, $0 (х) = $о(ж) = 0, а соответствующее S состоит из единственного элемента — множества {1, 2,..., г}. Приведем пример с г = 3, к = 2:
det М2 (у) = Р1 Р2Рз(Х2 - Ж1)2(жэ - Ж1)2(Ж3 - Х2)2,
det Мі(#, у) = д(Жі)^(Ж2)ріР2(Х2 - Хі)2+
+ д(жі)^(жз)рірз(хз - Хі)2 + д(ж2)д(жз)р2рэ(хз - Ж2)2
V* log det М2 (у) = [
Х1— Х2
+
Х1— Х2 Х2—Х1
+
Х2—хз
хз-хі
+
хз-Ж2
Р1
Р2
рз
где
Мх =
V2 logdet М2 (у) = Мх 0
0 Мр
2 2
(х1—Х2)2 (х1 хз)2 ' 2 (х1
(х1—х2)2 ‘ 2 (х2 2
(х1-хз)
МР = diag
(х1-хз) 1 '
(х2—хз)
р?:
р2:
рз
Конечные выражения для градиента и гессиана logdet М1(д,у) не приведены ввиду их размера. На практике вычислять их лучше не с помощью итоговых формул, а конструируя в соответствии с вышеприведенными утверждениями.
4. Обобщение результатов
Мы показали, что расчет Д,г в эквиваленте изложенного ранее алгоритма решения ЛМН-релаксации (5) может быть произведен без вычисления моментов у. Вследствие этого ограничения, связывающие порядки полиномов в задаче ПН с размером 66
2
2
2
2
т
2
2
2
2
2
2
1
1
неравенств и количеством неизвестных в ЛМН-релаксации, а также количеством атомов в новом алгоритме оптимизации, более не являются неотъемлемой частью метода решения.
Если трактовать не как [ 9 ], а как характеристику вза-
имодействия группы атомов с ограничением $ (х) ^ 0,7 мы сможем взять в качестве данных ограничений полиномы произвольного порядка и даже достаточно гладкие неполиномиальные функции. Аналогичное наблюдение относится и к целевой функции / (х). Таким образом, новый алгоритм естественным образом переносится на обобщение задачи (1) для более широкого класса функций. Конечно, гарантии нахождения глобального минимума для данной расширенной трактовки задачи и алгоритма оптимизации уже не будут иметь места: они по-прежнему будут относиться только к классу задач, для которых мы изначально строили ЛМН-релаксации. Тем не менее возможность одновременного исследования нескольких локальных минимумов в окрестности начального положения атомов — в том числе и минимумов, находящихся в других компонентах связности (см. раздел 5), — представляется достаточно ценной.
Расширенная трактовка задачи и алгоритма оптимизации лишает нас возможности рассматривать М*^($г,у) как функции моментов в соответствии с (3), так как количество коэффициентов д7, а следовательно, и необходимых моментов, теперь может быть неограниченным. Поэтому альтернативную высокоуровневую интерпретацию необходимо дать также и целевым функциям
т \
logdet М*(у) + ^ logdet М*—*($г,у) I , г=1 /
в которые входят данные матрицы. Рассмотрим компоненты данных функций по отдельности.
7В силу структуры det Mk-d^ (д, у), задаваемой утверждением 3, данная величина гарантированно останется положительной (а значит, будет выполняться и условие Мк-^ (д, у) > 0), если в процессе поиска не более ^ атомов одновременно окажутся на границе области {х|$;(ж) > 0}.
/(г)(у) = сТу -
• Слагаемое
СТУ = Е /[г]
г=0
г
г
Е^'
Е^
Ер- / (х-)
-=1
-=1
является взвешенным средним значением целевой функции по имеющейся группе атомов.
• В соответствии с утверждением 3,
Соответствующее слагаемое играет роль своего рода потенциальной функции, отталкивающей веса р от 0, а атомы друг от друга.
• Аналогичным образом logdet Мк_^(д*,у) играет роль потенциальной функции, дополнительно отталкивающей атомы от границы допустимой области |ж|дг(ж) ^ 0}. Отметим, что последнее относится к группе атомов в целом; индивидуальные же атомы вполне способны выходить в запрещенную зону.
Для того чтобы можно было вычислять Д,г с помощью формулы (10), по существу представляющей собой результат применения леммы 1 к функции с градиентом дх (дг) и модифицированным гессианом Нх (—), необходимо выполнение условий , ^ € со1 —; — ^ 0. Ранее эти соотношения следовали из € со1 и ^ 0 ив теореме 1 отдельно не доказывались. Поскольку новая трактовка алгоритма подразумевает отказ
ІС£ det Мд, (у) = 1о§ det Мй (1, у)
от использования вектора моментов, необходимо убедиться в выполнении данных условий без привлечения у и связанных с ним величин.
• В общем случае, если какая-либо из функций дДж) является полиномом степени выше 2di или неполиномиальной функцией, матрица Hz невырождена, и требование
, vz G col Hz выполняется. Отметим, что, несмотря на невырожденность, число обусловленности Hz может быть велико. Поэтому здесь целесообразно дополнительно применять метод Левенберга-Марквардта: использовать в (10) HHz + А/ вместо HHz, где / — единичная матрица; А — некоторый коэффициент.
• Следующая теорема, будучи применена к компонентам gz и Hz (11), доказывает неравенство Hz ^ 0, а также дает альтернативные формулы для вычисления gz и Hz, имеющие аналогичную (7) структуру. Особо отметим независимость модифицированного гессиана от вторых производных F(ж).
Теорема 3. Пусть F = F(ж) — симметричная вещественная матрица8 с дважды дифференцируемыми по ж элементами, и F(ж0) > 0. Пусть 9
= V logdet F(жо), Ях = Ях - ЯХ,
Яж = V2 logdet F (жо), Яу = ^ Fj т(жо)^^ (жо)).
ij
8Использование ж вместо z имеет целью соблюдение соответствия обозначений теоремам 1 и 2.
9 Несмотря на использование введенного ранее обозначения Яу, вектор у не фигурирует ни в формулировке, ни в доказательстве данной теоремы.
" Л
-
Доказательство. Далее, если не указано обратное, все функции и производные вычисляются в точке Ж0.
В силу симметричности Р найдется не зависящая от ж ортогональная матрица Q такая, что Р(ж) = ^ТС(ж)^ (при этом матрица С (ж) также симметрична), и С(ж0) — диагональная матрица (с положительными диагональными элементами). Тогда Ну = Хг • В самом деле,
№ >» = Е — ( аж^^) =tr (Р- ^Р) =
г,-
ажр ажд
Последнее выражение упрощается до ( ^г 0 йгг ) в силу диаго-
V Угг / рд
нальности С(ж0).
Представим определитель det Р как
det Р(ж) = /0(ж) - /1 (ж) + /2(ж),
/0 = Л дгг,
г
У1 ^ ' дг-д-г ^ ' ^г- /1,г-, /1,г- = Л ,
г<- &№,-} г<- &№,-}
/2 = Е (-1)*(а) П дгаг,
аЄА г
где A — множество перестановок, в которых более двух элементов находятся не на своих местах; N (а) — количество инверсий в а (/о соответствует тождественной перестановке в формуле Лейбница, /1 — перестановкам с одной инверсией, /2 — всем остальным). Поскольку gj (жо) = 0 при i = j, имеем:
• Определитель det F(жо) = /о(жо), так как /1(жо) = 0, /2(жо) = 0.
• Слагаемые, из которых состоят /1 и /2, содержат не менее двух (с учетом степени) внедиагональных элементов G в качестве сомножителей. Их производные, в соответствии с формулой для производной произведения
(uvw ...) = u'vw ■ ■ ■ + uVw ■ ■ ■ + uvw; ■ ■ ■ + . . . ,
содержат в каждом своем слагаемом не менее одного вне-диагонального элемента. Следовательно, V/^жо) = 0,
V/2(жо) = 0, и V det F(жо) = V/о^о).
• Для вторых производных, рассуждая аналогичным образом, получаем V2/^^ = 0.
Для градиента gx в точке жо получаем:
1 т-, V det F V/о „, „
gx = V logdet F = = —— = V log /о =
det F /о
= £ V log gii = £ V^,
дгг
г г
откуда
(gx )p = E # = * (G-1 (,1G)) =
= tr (qf-1 (diF) qt) =tr
Перейдем теперь к вычислению модифицированного гессиана. Найдем V2 log /о и V2 /1 в точке жо:
V2 log /о = £ V2 log gii = £ ( V2g“ <Vgii><Vg“>T
gi2i
ii
d2 2 t ( d \ f d \ ^ 2 d2^,
ix^xq 9j = 2 vdxp 9iv v"dxq 9iv ,ij +9j dxpdxq+
+2gij(( dxpTdxq 9ij) /Mj+
+ (gij) (dxq /Mj) + (dxqgij) (dx^ /Mj = 2 (d^9ij) (dxq 9ij) /mj = (2 (Vgij) (Vgij )T /Mj )Рд:
откуда
V2/l = 2 Е (V9ij)(V9ij)T /l,ij
i<j
= 2 £ (Vgij) (Vgy)T П gkk = 2/о £ (Vgt<Vgj)T.
i<j k0{ij} i<j гг
Собирая найденные элементы вместе, в точке жо получаем:
Ях = V2 log det F - ЯУ =
V2det F (Vdet F)(Vdet F)T £ V2gii-
det F (det F )2 ^ g..
V2/o - V2/l - (V/o)(V/o)T - е Vfe
o /o2 i gii
V2/o (V/o)(V/o)T v V2giiА _ V2/l
/o У,2 ^ 9гг / Л
V2/0 - Е ^-/ =
= ( ^ (Уугг)(Уугг)Т\ __ 2 ^ (Ууг.7) (Уугі)Т
угг / ,-<, угг
г<-
Т
г- ) ( У У г- )
(Ууг? ) (Ууг; )т
откуда ННХ ^ 0. Дальнейшие преобразования правой части приводят к финальной формуле:
(НУ ) = | ^ (Ууг-) (Ууг-)
(нж)рд = I
г-/ V V г-
х)рд — I / V
7- УггУ--
/ рд
угі) (~аідУгі) = ~ахру- ~аіду-
уггуіі угг —
Е (С-1 (С))- ((-С) ^
= - tr (С-1 ( -Ас) с-1 ( ,1с)) =
"жр ) V "жд ) /
= -QF-‘< "ГРР-‘ (£Р) ^) =
= - ^ ( Р-1 ( т^р) р-1 (
йжр / \ йжд
5. Пример
Рассмотрим задачу (1) с т = 1, /(ж) = ж2, у1(ж) = -2ж4 +
4ж2 - 1:
2
/ = Ш1П ж2,
2 -1 ^ 0.
- 2ж4 + 4ж2 - 1 ^ 0.
т
Графики / (ж) и ді(ж) приведены на рис. 1. Неравенство $і(ж) ^
0 задает область поиска в виде двух отрезков, на которых / (ж)
имеет по одному локальному минимуму в точках ±у^1 — ^2/2 ~ ±0,5412 соответственно.
2
1
1.5 -1.0 0.5 1.0 1.5
1 -2
Рис. 1. Графики / (ж) и д1(ж)
Решим данную задачу методом моментов. Размеры матриц в ЛМН-релаксации определяются параметрами ^ =
|" 1 degд1(ж)] =2 и к ^ max{deg/(ж),^1} = 2. Согласно [7], мы должны решить ЛМН-релаксацию для двух последовательных значений к и, убедившись, что ранг получаемых матриц моментов совпадает, извлечь координаты минимумов согласно изложенному в данной работе алгоритму. В рассматриваемом случае достаточно взять к = 2 и к = 3.В частности, для к = 3 соответствующая ЛМН-релаксация (5) имеет вид
/ * = пут у[2],
Мз(у) =
у[0] У[1] У[2] у[з]
У[1] У[2] у[з] У[4]
У[2] у[з] У[4] У[5]
- у[з] У[4] У[5] У[6]
_2У[4] + 4У[2] —
_2У[5] + 4У[з] — У[1]
—2У[5] + 4У[з] — У[1]
— 2У[6] + 4У[4] — У [2]
= 1,
и решение у* = [У[0] ... У[6]] =[1 0 0,2929 0 0,0858 0 0,0251]. Ранг Мз(у*) равен 2, и извлеченные из нее координаты двух экстремумов совпадают с указанными выше.
74
Продемонстрируем теперь процесс решения той же задачи предлагаемым в данной работе методом. Как указано в разделе 4, здесь у нас больше свободы в выборе к и di. Рассмотрим несколько вариантов.
1) Пусть к = di = 0, а количество атомов r = к + 1 = 1. Тогда
z = [xi pi]T,
Mo(y) = [ pi ] , Mo(gi,y) = [ pi(-2x1 + 4xf - 1) ] ,
и целевые функции f(i) имеют вид
f w(z) = pixi - ^(i) (logpi + log (pi(-2xt + 4xf - 1))) .
В зависимости от того, какому из составляющих область поиска отрезков принадлежит начальное приближение, мы найдем один из двух экстремумов задачи. На рис. 2 показана траектория xi для начального приближения xi0) = 1, P(i0) = 1.
l.O
O.5
-O.5
-l.O
1O 2O 3O 4O!
Рис. 2. График ж1 для к = 0, ^1 =0. Серым цветом отмечена область у1 (ж) ^ 0
X
Здесь и далее алгоритм делает серию из 15 шагов поиска с = 1, а после нее — 5 серий по 5 шагов, уменьшая ^ в 4 раза в каждой новой серии:
/°"Л4> = 1, ^(15"Л9) = 1/4, ..., ^з5...з9) = 1/44.
(Из графика очевидно, что для выбранного значения к такое количество итераций является излишним; оно использовано здесь исключительно для единообразия масштаба.)
Таким образом, действие предложенного алгоритма в данном случае эквивалентно простому локальному поиску с использованием логарифмических барьерных функций и ньютоновским выбором направления — за исключением того, что мы используем модифицированный гессиан Н, который всегда положительно полуопределен.
2) Пусть к = ^1 = 1, г = 2. Тогда г = [ж1 ж2 р1 р2]Т,
Р1 + Р2 Р1ж1 + Р2ж2
22
PlXl + Р2Х2 PlX2 + Р2Х2
Ml(y) =
Mo(gl,y) = [ Pl(-2x4 + 4x1 - l) + P2(-2^4 + 4x2 - l) ]
f(i) (z) = PlX2 + P2x2 - (log (P1P2(X2 - Xl)2) +
+ log (pl(-2x4 + 4x2 - l) + P2(-2x2 + 4x2 - l))) .
Выберем в качестве начального приближения случайный вектор z(0) такой, что xl0) и x20) находятся вблизи точки
l, а веса принадлежат интервалу (О; l) и удовлетворяют ограничению pl_0) + p20) = l. Траектории xl и x2 показаны на рис. За; pl и p2 — на рис. Зб. Особого внимания здесь заслуживает своеобразный «туннельный эффект»: один из атомов переходит через запрещенную область (gl(x) < 0) в соседний допустимый отрезок, за счет чего мы находим оба экстремума. Последнее становится возможным при достаточно существенном уменьшении веса атома, позволяющем матрице Mo(gl,y) оставаться положительно определенной.
3) Дальнейшее увеличение количества атомов приводит к тому, что они разделяются на две группы:
• Атомы, попавшие в глобальные экстремумы, приобретают веса на интервале (0; l). При этом возможно попадание нескольких атомов в окрестность одного экстремума.
и
а б
Рис. 3. Графики Жг и р для к = 1, = 1
• Веса прочих атомов, в том числе попавших в неглобальные экстремумы (отсутствующие в данной задаче) или вышедших за пределы допустимой области, стремятся к 0 при уменьшении ^(г).
На рис. 4 показан пример траекторий жг и рг для к = 2,
^1 = 1.
аб Рис. 4. Графики жг и рг для к = 2, ^1 = 1
X
При росте количества атомов растет объем вычислений, а также — как видно из графиков — замедляется сходимость алгоритма. Поэтому, как и в методе моментов, в новом алгоритме целесообразно использовать как можно меньшие значения к. Несмотря на это сходство, собственно критерии выбора к отличаются: если в методе моментов мы исходили из порядков поли-
77
номов f (ж) и gi(x), то в новом алгоритме основным фактором является ожидаемое количество экстремумов в области поиска.
6. Заключение
В данной работе представлена трансформация основанного на теории моментов метода решения оптимизационных задач с полиномиальными целевой функцией и ограничениями. Данная трансформация возвращает модифицированную задачу оптимизации («ЛМН-релаксацию») из пространства моментов обратно в пространство поиска исходной задачи, расширенное путем рассмотрения нескольких взаимодействующих «атомов» и их весов вместо единственной точки-аппроксимации. Это позволяет устранить наиболее существенные недостатки метода ЛМН-релаксаций, а также расширить круг решаемых задач.
В статье детально рассмотрено применение нового метода к задачам одномерной оптимизации. Последующие работы будут посвящены распространению полученных результатов на задачи многомерной оптимизации и системы полиномиальных матричных неравенств.
Литература
1. APKARIAN P., TUAN H.D. Parameterized LMIs in control theory // SIAM J. on Control and Optimization. - 2000. -Vol. 38, №4. - P. 1241-1264.
2. BLONDEL V., TSITSIKLIS J. NP-hardness of some linear control design problems // SIAM J. on Control and Optimization. - 1997. - Vol. 35, №6. - P. 2118-2127.
3. BOYD S., GHAOUI L.E., FERON E., BALAKRISHNAN V. Linear matrix inequalities in system and control theory. -SIAM, 1994.
4. BOYD S., VANDENBERGHE L. Convex optimization. -Cambridge University Press, 2004.
5. CURTO R.E., FIALKOW L.A. The truncated complex K-moment problem // Trans. AMS. - 2000. - Vol. 352, №6. -P. 2825-2856.
6. DELIBASI A., HENRION D. Hermite matrix in Lagrange basis for scaling static output feedback polynomial matrix inequalities // International J. of Control. - 2010. - Vol. 83, №12. - P. 2494-2505.
7. HENRION D., LASSERRE J.-B. Detecting global optimality and extracting solutions in GloptiPoly // Positive polynomials in control. - 2005. - P. 1-18.
8. HENRION D., LASSERRE J.-B. Convergent relaxations of polynomial matrix inequalities and static output feedback // IEEE Trans. Automatic Control. - 2006. - Vol. 51, №2. -P. 192-202.
9. HENRION D., LASSERRE J.-B., LOFBERG J. GloptiPoly 3: moments, optimization and semidefinite programming // Optimization Methods and Software. - 2009. - Vol. 24, №4-
5. - P. 761-779.
10. HENRION D., LOFBERG J., KOCVARA M., STINGL M. Solving polynomial static output feedback problems with PENBMI // Proc. of the 44th IEEE Conference on Decision and Control. - 2005. - P. 7581-7586.
11. IWASAKI T., SKELTON R. Parametrization of all stabilizing controllers via quadratic Lyapunov functions // J. of Optimization Theory and Applications. - 1995. - Vol. 85. -P. 291-307.
12. http://www.penopt.com (дата обращения: 01.11.2011).
13. KOCVARA M., STINGL M. PENNON: A code for convex nonlinear and semidefinite programming // Optimization Methods and Software. - 2003. - Vol. 18, №3. - P. 317-333.
14. LASSERRE J.-B. Global optimization with polynomials and the problem of moments // SIAM J. on Optimization. - 2001. -Vol. 11, №3. - P. 796-817.
15. PETERSEN K.B., PEDERSEN M.S. The
matrix cookbook. — 2008. - URL:
http://citeseer.ist.psu.edu/viewdoc/summary?doi=10.1.1.113.6244 (дата обращения: 01.11.2011).
ATOMIC OPTIMIZATION, PART 1: SEARCH SPACE TRANSFORMATION AND ONE-DIMENSIONAL PROBLEMS
Vladimir Pozdyayev, Arzamas Polytechnical Institute of R. E. Alekseev Nizhny Novgorod State Technical University,
Arzamas, Cand.Sc., assistant professor (vpozdyayev@gmail.ru).
Abstract: Optimization problems with polynomial objective function and inequality constraints are considered. A transformation of the moment theory-based solution method is presented, which allows to construct equivalent solution algorithms working in the augmented original search space instead of the moment space. One-dimensional problems are analyzed in detail.
Keywords: nonlinear programming, matrix inequalities, polynomial inequalities, moment theory.
Статья представлена к публикации членом редакционной коллегии П. С. Щербаковым
