CC BY
68
УДК 681.5.011 +519.853.4 ББК 32.965 + 22.18
АТОМНАЯ ОПТИМИЗАЦИЯ, ЧАСТЬ 2: МНОГОМЕРНЫЕ ЗАДАЧИ И ПОЛИНОМИАЛЬНЫЕ МАТРИЧНЫЕ НЕРАВЕНСТВА1
Поздяев В. В.2
(Арзамасский политехнический институт (филиал) Нижегородского государственного технического университета им. Р. Е. Алексеева, Арзамас)
Рассмотрены многомерные задачи оптимизации с полиномиальной целевой функцией и ограничениями в виде полиномиальных матричных неравенств. Представлена трансформация основанного на теории моментов метода их решения, позволяющая существенно снизить его вычислительную сложность, сохранив способность решать задачи интересующего нас класса.
Ключевые слова: нелинейное программирование, матричные неравенства, полиномиальные неравенства, теория моментов.
Введение
Рассмотрим задачу нахождения глобальных экстремумов полиномиальной целевой функции на множестве, заданном полиномиальными неравенствами (ПН):
/ * = шш / (х),
X
(1) §г(х) ^ 0,
х е Мга, г = 1,..., т,
1 Работа выполнена при финансовой поддержке РФФИ, грант №1208-31440.
2 Владимир Васильевич Поздяев, кандидат физико-математических наук, доцент (vpozdyayev@gmail.com).
или полиномиальными матричными неравенствами (ПМН):
/ * = шт / (х),
(2) Сг(х) > 0,
х е Мга, г = 1,..., т, где /(х) и §г(х) — (не обязательно выпуклые) полиномы; О^х) — матрицы, элементы которых являются полиномами от х, а знак неравенства в (2) понимается как требование положительной по-луопределенности. Далее данные задачи мы будем называть соответственно задачами ПН и ПМН.
В первой части данной статьи [1] рассматривался предназначенный для решения таких задач метод глобальной оптимизации [4, 5, 7], фундамент которого образуют теория разложения полиномов в сумму квадратов и двойственная ей теория моментов. Данный метод позволяет найти все глобальные экстремумы, используя для этого сведение исходной задачи к иерархии систем линейных матричных неравенств (так называемых ЛМН-релаксаций), решение которых представляет существенно меньшую трудность. Тем не менее, применимость данного метода ограничена двумя факторами:
• комбинаторным взрывом количества неизвестных и размера матриц в ЛМН-релаксациях;
• ухудшением обусловленности составляющих ЛМН-релаксации матриц и последующим снижением точности расчетов по мере рассмотрения последовательности ЛМН-релаксаций.
Обе проблемы вызваны тем, что в ЛМН-релаксациях пространством поиска является пространство не переменных исходной задачи, а их моментов соответствующих порядков. Для уменьшения влияния данных проблем на результат авторы метода применяют вспомогательные техники, такие как представление исходной задачи в форме с как можно меньшим количеством неизвестных; уменьшение количества моментов-переменных в ЛМН-релаксациях за счет использования линейных связей между ними
(если постановка задачи содержит зависимости такого рода); масштабирование переменных таким образом, чтобы искомые экстремумы удовлетворяли условию ||х*|| ^ 1.
Существуют иные алгоритмы, основанные на подходе аналогичного вида (построение иерархии аппроксимаций исходной задачи и решение их более или менее стандартными методами). Возникающие при этом структуры в общем случае адаптированы к классу решаемых задач и могут иметь сниженную вычислительную сложность за счет усиления консерватизма аппроксимаций; см., например, сравнение методов решения задач бинарного программирования в [9]. Кроме того, модификация отдельных этапов метода [7] (в первую очередь методики построения аппроксимаций) может позволить решать задачи самого разного вида, в том числе имеющего неочевидное отношение к задачам ПМН (общая структура метода при этом остается неизменной); см., например, [3, 6, 8].
В [1] был предложен альтернативный подход, основанный на исключении из процедуры решения пространства моментов. Это достигается за счет дальнейшей трансформации ЛМН-релаксаций с целью возвращения задачи в (расширенное) исходное пространство поиска. Также была разработана вычислительная схема, позволяющая с минимальными изменениями применить к новой задаче метод внутренней точки в прямой форме, который может использоваться для решения ЛМН-релаксаций. Данные вопросы детально рассматривались для одномерных задач оптимизации в форме ПН.
Данная статья посвящена дальнейшему развитию техники преобразования пространства поиска применительно к многомерным задачам оптимизации и полиномиальным матричным неравенствам. Раздел 1 содержит сведения о базовом методе, относящиеся к такого рода задачам. В разделе 2 приведены результаты из [1], необходимые для дальнейших построений. Основной раздел 3 распространяет данные результаты на задачи ПМН вида, характерного для теории управления. В разделе 4 приведены примеры применения полученных результатов.
1. Базовый метод
Основные положения метода решения задач ПН (1), опубликованного в [7], были изложены в предыдущей части статьи. Приведем необходимые нам в дальнейшем элементы.
Пусть Ьг(х), х Є К”, — вектор, состоящий из одночленов, образующих базис пространства многочленов порядка не выше г:
Ьг (х) = [ 1 жі х2 ... хп х1 хіх2 ...
... хіхп х2хз ... х” ... хі ... х” ]Т ,
а зп(г) = СП+Г = (п+гГ — его размерность. Каждому одночлену из Ьг (х) поставим в соответствие вектор а Є М”, ^і а ^ г (далее будем записывать как а ^ г), показателей степеней х1, х2, ..., хп; обозначим ха = х^1 х^2 ...х^". Для произвольного вектора р Є К5п(г), ассоциированного с пространством моментов х степени не выше г, будем индексировать его элементы двумя взаимозаменяемыми способами: по номеру элемента и по вектору показателей степеней; порядок элементов будем считать соответствующим структуре Ьг (х). Таким образом, Р = Ык^п(г) = [Р«]«ЄН", а<г, в том числе рі = Р[0,0,...,0], Р2 = Р[1,0,...,0] и т. д. Аналогичным образом будем индексировать строки и столбцы матриц там, где это применимо.
Рассмотрим некоторую (неизвестную) меру ^ и соответствующий ей вектор моментов у:
у = J Ьг (х) ф.
Пусть б,і = |"І deg^і(х)], а к удовлетворяет ограничениям 2к ^ deg /(х), к ^ ^. Пусть /а — коэффициенты /(х) в базисе Ь2к(х), так что
[ /(х^ = [ ^ /аха d^ = ^ /аУа.
^ а^2к а^2к
ЛМН-релаксацией (1) будем называть задачу
/* = ШІП V /аУа,
У '
У
а^2к
(3)
Мк(у) ^ 0,
Ми-^(&,у) ^ о, г = 1,..., т, у[0,0,...,0] = 1
где матрица моментов М^(у) и локализующие матрицы Ык-а1 ($г,у) конструируются исходя из соотношений
В [7] (теорема 4.2) показано, что, с учетом некоторых непринципиальных ограничений, при к ^ то величина экстремума ЛМН-релаксации стремится к величине экстремума исходной задачи ПН. Более того, как правило, уже при конечных (и относительно небольших) значениях к данные величины становятся равны, а вектор моментов решения задачи ПН является решением соответствующей ЛМН-релаксации. Достаточным условием достижения такого значения к является
где у* — решение ЛМН-релаксации, а й = шах^ ^. Если оно выполняется, то у* представляет собой вектор моментов г-атомной меры3, атомы которой х*7, ^ = 1,..., г, соответствуют глобальным минимумам (1). Данные атомы могут быть извлечены из у* путем решения системы полиномиальных уравнений (в которую для г-атомных мер превращается (4)) с помощью алгоритма, представленного в [4].
3 N-атомная мера — мера, носитель которой является множеством
из N точек (атомов).
(4)
(5) Мк-^(#,у)= / Ьк-й(х)Ьк-й(х)Т#(х^.
г = гапк Мк (у*) = гапк Мк-^(у*),
Одним из способов решения задач ПМН (2) является преобразование их к форме ПН с последующим применением описанного метода. Но более эффективным подходом является построение ЛМН-релаксаций непосредственно для исходной матричной формы неравенств способом, описанным в [5]. А именно: вместо (5) будем конструировать локализующие матрицы, исходя из следующего соотношения:
Мй-й(С, у) = J (Ь^(ж)Ь^(ж)т) ® С(ж) ф. Например, для п = 2, к — d = 1 и
С(ж) =
локализующая матрица имеет вид
ж1 2
2 Ж2
М2(С,у) = У А ф = В,
А =
В =
1 2 іЖ 2Ж 2 ю ж2 2ж1 2ж1 ж1ж2 ж1ж2 2ж2 2ж2 ж22
ж2 2ж1 2ж1 ж1ж2 /у»3 О'У»2 ж 1 ^ж 1 2ж1 ж1ж2 3^01 О! ж1 ж ж1 ж1 2ж 22 ж2 ж12 21 1 ж ж2
ж1 ж2 2ж2 _ 2ж2 ж2 У[1,о] 2У[о,о] 2У[о,о] У[о,1] ж2ж2 2ж1ж2 2ж1ж2 ж1ж2 У[2,о] 2У[1,о] 2У[1,о] У[1,1] /у» /у»2 О'У’2 ж1 ж2 ж2 О/у»2 /у»3 ж2 ж2 У[1,1] 2У[о,1] " 2У[о,1] У[о,2]
У[2,0] 2У[1,о] 2У[1,о] У[1,1] У[з,о] 2У[2,о] 2У[2,о] У[2,1] У[2,1] 2У[1,1] 2У[1,1] У[1,2]
У[1,1] 2У[о,1] - 2У[о,1] У[о,2] У[2,1] 2У[1,1] 2У[1,1] У[1,2] У[1,2] 2У[о,2] 2У[о,2] У[о,з] .
2. Результаты части 1
Результаты, полученные в предыдущей части, были посвящены двум вопросам: определению эквивалентного направления 100
поиска в методе внутренней точки при трансформации пространства поиска; применению полученного результата к одномерным задачам ПН. Приведем подробности в форме, наиболее подходящей для дальнейших построений.
2.1. ЭКВИВАЛЕНТНОЕ НАПРАВЛЕНИЕ ПОИСКА Пусть задача ПМН
/ * = шш / (х),
X
(6) ^Кх) ^ o,
рХх = 1,
х е Мга, г = 1,..., т,
может быть получена из задачи ЛМН
/*& т
/ = шт с у,
у
^г(у) ^ 0,
гр
^уУ = 1
(7)
у
у е Мга, г = 1,..., т, с помощью замены4 (трансформации пространства поиска) у =
у(х).
Далее мы будем рассматривать алгоритм решения (7) на основе метода внутренней точки в прямой форме с ньютоновским направлением поиска, представляющий собой серию подзадач минимизации целевых функций вида
т
/(г)(у) = сТу — ^ logdet ^ (у^
.7=1
где |^(г)} — монотонно невозрастающая сходящаяся к 0 вещественная последовательность. Каждая подзадача решается на
4 Одним из условий результатов [1] было отличие от 0 якобиана у(х) внутри допустимой области. Здесь мы опускаем это условие, поскольку финальные результаты целенаправленно конструируются так, чтобы не содержать упоминаний о данной трансформации.
ограниченном неравенствами р(у) ^ 0 отрезке прямой, проходящей через у(г) и имеющей предписываемое методом Ньютона направление, которое с учетом ограничения в виде равенства имеет вид
(8) ду«у(,))=н- (—ду+Цтт ">) ■
где н- — произвольная5 обобщенная обратная к Ну матрица;
т
ду = Vy/(%(г)) = с — £ Уу logdet р(у«),
.7 = 1
Ну = Уу/(%«) = -^ £logdet р(/)),
.7 = 1
а элементы слагаемых под знаками сумм могут быть найдены как
(уу logdet р(у))г = ^ (Р-1(у) (^Р(у))) ,
(Уу logdet р(у))г7 = — ^ (Р-1(у) (-дутР(у))
^-1 (у) (^ ^ (у))) .
Как показано в предыдущей части, такой алгоритм естественным образом переносится в пространство поиска задачи (6). При этом вспомогательные целевые функции имеют вид
т
(9) /(г)(х) = /(х) — ^(г) £ logdet р (х),
7=1
а направление поиска, эквивалентное в малом направлению (8), равно (см. теоремы 1, 2, 3 в [1] и комментарии к ним)
(10) Дх(х«) = Н- ( —дх + .
\ ^Т НX ^х /
5 Свободно выбирается из множества обобщенных обратных матриц при первом упоминании. Все дальнейшие вхождения этой матрицы имеют то же значение.
Здесь Н— — произвольная обобщенная обратная к Нх матрица, а градиент дх и модифицированный гессиан Нх находятся по формулам
т
дх = Ух/(%«) = Ух/(х«) — £ Ух logdet Р,(х(г)),
7=1
т
Нх = Vх/(%«) = —^« £Ух logdet Р,(х«),
7=1
где Vх = Ух — ^г(Уху»)"дуг; слагаемые под знаками сумм могут быть вычислены следующим обр(азом: ( ))
(Ух logdet Р(х))* = tr (V-1 (х) (-дхР(х))) ,
(11) (ух logdet Р(х))*7 = — ^ (Р-1(х) ( ахР(х))
Р-1(х) ( Ж-Р(х))) .
Мы видим, что (10) не зависит от вида трансформации у = у(х) и компонентов задачи (7). Более того, даже если не существует линеаризующей трансформации указанного вида6, полученный результат сохраняет свою применимость, хотя и теряет теоретические гарантии нахождения глобального минимума итоговым алгоритмом. Особо отметим сохранение при этом важного для метода Ньютона свойства Нх ^ 0, а также возможность использования в трансформации задачи (3) значений к и di, не ограничиваемых более порядками полиномов исходной задачи ПН/ПМН.
2.2. ОДНОМЕРНЫЕ ЗАДАЧИ
Указанные результаты были применены к задачам ПН с п =
1. В качестве трансформации пространства поиска был взят переход из пространства атомов в пространство моментов у = у(^).
6 Как видно из раздела 1, для задач ПН/ПМНлинеаризующая трансформация найдется всегда, но минимально допустимая размерность у при этом может быть больше п. Если же в задаче фигурируют неполиномиальные функции, такой трансформации с вектором у конечной размерности может и не быть.
Здесь г — вектор в пространстве атомов МГх(га+1) = М2г, имеющий структуру
г = [ Ж1 Х2 ... Жг Р1 Р2 ... Рг ]Т ,
где Жі Є М — атомы; рі є (0; 1) — их веса; вектор моментов у имеет размерность ^(2к) + 1 = 2к + 2 и состоит из элементов вида
г
У, = У[і-і] = £РіЖі—1, 2 = 1 2,..., 2к + 2.
і=1
Старший, (2к + 1)-й момент не входит ни в (3), ни в Ду, и нужен исключительно для равенства размерностей у и г, которое обеспечивается дополнительным соотношением г = к + 1.
Интерпретируя ЛМН-релаксацию рассматриваемой задачи как форму (7), мы можем записать эквивалентную задачу вида (6) без использования вектора моментов следующим образом:
(12)
Г
/ * = тіп £ р, / (жі ^
* і=1 Г
ВД = £Рі (^(жі^(жі)Т) ^0, і=1
Г
Я (г) = £ р, (^-^ (ж, )Ьк-гіі (ж, )Т) ді(жі) ^ 0, і = 1,...,т, і=1
^Тг=£ р,=1.
І=1
Данная форма задачи напрямую подходит для решения модифицированным методом внутренней точки с вычислением эквивалентного ньютоновского направления поиска Дг(г(і)) по формуле (10).
3. Многомерные задачи
Задачи ПН/ПМН в общем случае являются NP-трудными [2], и предлагаемый способ их трансформации данную проблему не решает. С другой стороны, благодаря смене пространства поиска, вычислительная сложность процедуры решения задачи определяется теперь не столько формальными характеристиками последней, сколько фактической сложностью структуры области поиска и целевой функции. Например, даже если задача заведомо является выпуклой, но имеет относительно большое количество неизвестных и задействует полиномы не слишком малого порядка, оригинальный метод вынужденно столкнется с вышеупомянутым комбинаторным взрывом размера ЛМН-релаксаций. С помощью же новой схемы мы сможем сконструировать алгоритм с единственным атомом, реализующий локальную оптимизацию и требующий существенно меньших вычислительных ресурсов.
Далее мы будем ориентироваться на задачи родом из теории управления, зачастую изначально представляемые в виде ПМН и имеющие невысокий порядок полиномов, не слишком малое количество неизвестных, а также область поиска с относительно несложным характером невыпуклости (см. пример в разделе 4.2)7. Последнее условие подразумевает, что невыпуклость как таковая не приводит к наличию большого количества ложных локальных минимумов, сложному рельефу вспомогательных целевых функций и т. п. Таким образом, мы можем рассматривать невыпуклость задачи не как бинарный фактор, вынуждающий нас при его наличии в лучшем случае искать альтернативные формулировки задачи с более приемлемыми количественными характеристиками, а в худшем — смиряться с катастрофическим ростом объема вычислений — как неявную характеристику,
7 Примером задачи с патологическим характером невыпуклости, который мы не рассматриваем, может служить любая задача оптимизации на множестве битовых векторов (такие задачи представимы в виде ПН, поскольку условие х^ € {0; 1} эквивалентно системе х^ > 0; х4 < 1; хДх* — 1) > 0). В этом случае область поиска предсталяет собой дискретное множество из 2П точек.
под которую подбираются параметры алгоритма поиска решения (в первую очередь к и г).
Построим алгоритм поиска экстремума, основанный на описанной выше схеме и удовлетворяющий следующим требованиям:
1) он должен быть совместим как с задачами ПН, так и с задачами ПМН;
2) он должен избегать эффекта комбинаторного взрыва;
3) он должен позволять более тонко контролировать объем вычислений, в частности, выбирать необходимое количество атомов в зависимости от ожидаемого характера невыпук-лости задачи и количества локальных экстремумов.
Комбинаторный взрыв является неизбежным результатом повышения порядка ЛМН-релаксации в оригинальном методе, и попытка конструирования эквивалентной трансформированной задачи дала бы в этом случае аналогичный эффект. Поэтому дальнейшие построения основаны на ЛМН-релаксациях с минимальными нетривиальными значениями к и ^. Такие конфигурации, соответствующие первым двум требованиям, рассмотрены в разделе 3.1.
В разделе 3.2 представлено обобщение вычислительной схемы, позволяющее дополнительно уменьшать количество атомов до любого нужного значения, включая 1 (в соответствии с требованием 3).
Раздел 3.3 описывает способ компенсации консерватизма вычислительной схемы, вызванного отказом от использования ЛМН-релаксаций высоких порядков.
3.1. БАЗОВАЯ КОНФИГУРАЦИЯ
По аналогии с одномерными задачами рассмотрим пространство атомов, состоящее из векторов вида
где Жг € К” — атомы, рг £ (0; 1) — их веса. Соответствующий вектор моментов имеет вид у = у(г) = ^Г=1 ргЬд,(Жг).
Подставив у = у(г) в релаксацию (3) задачи (2), получаем аналог (12), отличающийся структурой векторов г (имеющего указанный выше вид), уг (который теперь состоит из гп нулей и г единиц), Ьд(ж), (ж), а также матриц Рг(г):
Базовая версия предлагаемого далее подхода использует минимальное нетривиальное значение к = 1 в сочетании с минимально возможным количеством атомов, дающим невырожденную матрицу Ро(-г): г = зп(к) = п + 1. Величины I при этом могут принимать значения 0 и 1. Определим, какие из данных значений целесообразно использовать. Отметим, что Ро(^) можно формально считать разновидностью РДя) с 10 = 0 и д0(ж) = 1.
Утверждение 1. Пусть к = 1, г = 5п(к) = п + 1, и С(ж) є М1х1. Определитель матрицы
(13)
Г
Р(я) = Р' (6к-^і(жІ)Ь^і(ж'Г) ® ^і(жІ).
'=1
Г
Р(я) = £Рі (Ьк-гі(жі)Ьк-гі(жг)Т) ® С(жі)
равен
при I = 0 и
Г
det Р (я) =ёе^£ рг С(жг)
при й = 1. Здесь V € К(га+1)хг —матрица, г-й столбец которой равен Ь1(жг) (п-мерная матрица Вандермонда порядка 1 для векторов жг):
1 1 ... 1
жц ж2і . . . жг1
V = жі2 ж22 . . . жГ2
ж1га ж2га . . . жгп
Доказательство. Формула для й = 1 следует из Ь0(ж) =
[1]. Пусть теперь й = 0; Н = Нт — произвольный квадратный корень из РгС(Жг); W = = diag(Hl, Я2,..., Н); / € К1х1 —
единичная матрица. Тогда
Г
Р(г) = £ Рг (Ьй-й(Жг)Ьй-й(Жг)т) ® С(ж*) =
г=1
Г
= £ (Ьй-й(Жг) ® Нг) (Ьй-й(Жг) ® Нг)Т = г=1
= ((V ® /^)((V ® /)W)т =
= (V ® /г)(WWT)(V ® /г)т,
так что (поскольку V ® /г и WWт — квадратные матрицы)
det Р(г) = det(WWT)(det(V ® /г))2 =
= (detdiag (р1С(ж1),р2С(ж2),... ,рГС(жг))) (det V)2г,
откуда получаем доказываемую формулу.
Мы видим, что при использовании йг = 0, г = 1, 2,..., т, необходимым условием сохранения положительной определенности Рг(г) в процессе решения задачи методом внутренней точки является положительная определенность каждой из матриц Сг(ж^) по отдельности. Данное требование лишает атомы ж^- возможности покидать допустимую область — что является одним 108
из ключевых аспектов алгоритма, особенно при исследовании областей поиска с несколькими компонентами связности. По этой причине использование таких значений ^ представляется нецелесообразным. Далее мы сохраняем ^о = 0, но полагаем все остальные ^ равными 1. Задача (6) тогда приобретает вид
Г
/ * = тт £ р- / ^
* -=1
Ро(-г) = Vdiag(pl,p2,... ,Рг)VT ^ 0,
(15) Г
^(*0 = £Р-^г(ж-) ^ 0, І = 1, . . . , Ш,
-=1
^ = £ Р- = 1.
-=1
3.2. РЕДУКЦИЯ
Представленная в предыдущем разделе схема имеет жесткое ограничение на количество атомов: г = (1) = п + 1. Между
тем, в реальных задачах ПМН данное число может быть довольно велико. Покажем, что выведенные ранее формулы допускают обобщение на случай г < п + 1.
Единственным элементом (15), несовместимым с изложенной схемой при г < п + 1, является неравенство Ро(^) ^ 0 (поскольку матрица Р0(г) в этом случае вырождена). Возвращаясь к интерпретации logdet Ро(^) как потенциального поля, отталкивающего атомы друг от друга (см. [1], раздел 4), отметим, что при г = 1 данное неравенство теряет смысл и его можно исключить. Пусть теперь г > 1. Определим воздействие на атомы аналогичного поля, действующего в проходящей через ж1, ж2, . . . , жг (г — 1)-мерной гиперплоскости X. Далее будем предполагать, что конфигурация атомов жь ж2, ..., жг не является вырожденной: они не лежат на гиперплоскости меньшей размерности.
Введем на указанной гиперплоскости ортонормированный базис, и пусть Мо Є Мгах(г-1) — матрица, столбцы которой являются элементами данного базиса. Эта матрица задает линейный
оператор, ставящий векторам х Є X в соответствие их координаты в новом базисе: X = Мтх (без ограничения общности будем считать началом координат в X решение системы Мтх = 0, х Є X). Пусть Хі = Мт Хі, і = 1,2,..., г. Обозначим также Ш = diag(p1,p2,... ,Рг), так что Ро(^) из (15) будет иметь вид Ро(-г) = УШУт.
Примем в качестве эквивалента Ро(^) на X матрицу
*Ъ(г) = V ШУт = М^тУШУтМ^ = М^т Ро(<г)М^,
где
V =
1
Х11
1
Х21
х1,г-1 х2,г- 1
1
хг1
Хг г-1
М^ = ^([1],М0) е М(га+1)хг.
Тогда
det ^0(2) = det(M0TV) ) det(VTM0) =
= ) det(VТМ0М^) = ёе^УТР'V),
где Р' = diag([1],P); Р = М0МТ — матрица проекции на пространство столбцов М0, которую можно найти, построив М0, или же как Р = М(МТМ)-1Мт, где столбцы М е Мгах(г-1) образуют произвольный базис X: например, х2—х1, —х1,..., хг — х1.
Отметим, что матрицы V, М0ТУ и Р'V имеют полный ранг по столбцам.
Выражение det(VТР'V) можно упростить. Для г = 1 данный определитель равен 1. Если г > 1, пусть
К
1 1 1 1 1 . -1 ■
0 1 0 .. .0
0 0 1 .. .0
0 0 0 .. .1
Є Мг
Тогда
) =ёе^Ят)ёе^У ТР/ТР'У ^(К) =
= det ((Р^К^Р^К)) .
Несложно видеть, что Рж1 линейно зависит от Р (ж2 — Ж1) = ж2 — Ж1, ..., Р (жг — ж1) = жг — ж1, а следовательно, существует такая нижняя треугольная матрица Ь с единичной диагональю (описывающая вычитание из первого столбца Р'VК остальных с соответствующими коэффициентами), что, домножив на нее справа матрицу
Р 'УК =
1 0 0 '
Р 1 ж2 — ж1 ж1 — ж
мы получим
Р 'УКЬ = ё1аё([1],М),
М = [ ж2 — ж1 | жз — ж1 | ... | жг — ж1 ] ,
так что ёе)((Р'УК)Т(Р'УК)) = det ((Р'УКЬ)Т(Р'УКЬ)) = det (МТМ). Финальная форма является грамианом системы векторов ж2 — ж1, ж3 — ж1,..., жг — ж1, и, таким образом, det(VTP'V) представляет собой квадрат (г — 1)-мерного объема параллелотопа, построенного на данных образующих (что согласуется с (14), поскольку при г = п + 1 определитель det V равен ориентированному объему параллелотопа с образующими того же вида). Таким образом,
(16) det —0(г) = det (МТМ) .
Найдем теперь компоненты logdet -^(г) и
V2 logdet -0(г). Будем считать X, М0 и М константами (так что 1 М0(г) = 0 и 1М(г) = 0): это позволит избежать зависимости модифицированного гессиана V2 logdet Р0(.г) от выбора базиса X (напомним, что данный гессиан определяется не только инвариантной по отношению к базису функцией logdet .-0(г), но и видом самой матрицы -0(г)). Мы можем
111
это сделать, поскольку нашей целью на данном этапе является определение воздействия барьерной функции, индуцированной неравенством —0(2) ^ 0 на X, на конкретную (текущую)
конфигурацию атомов.
Теорема 1. Пусть Р', V, Ш, —0(г) и —0(2:) —матрицы указанного выше вида, и = (Р^)- — произвольная обобщенная обратная к Р'V матрица, и С(г) = Р'иТШ-1иРТогда
(V logdet —0(2))* = tr (5(2) (—0(2))) ;
(V2 logdet —0(2))у = — tr (ОД (^5--0(г)) ОД (-0(г))) .
Доказательство. Пусть М0 — матрица указанного выше вида. Поскольку —0(2) — невырожденная матрица, существует единственная обратная к ней матрица (5(2), которая может быть найдена как 5(2) = МТиТШ-1иМ0. Это можно видеть, подставив данное выражение, а также —0(2) = M0TVWVТМ0, в —0(2)5(2)—0(2) и упростив результат с учетом равенств М0Л^(0Т = Р' = Р'Т и (P'V)-(Р'V) = I (в силу пол-норанговости Р^ по столбцам). Результатом данных действий является выражение, идентичное —0(2), а следовательно, (5(2) является обобщенной обратной — и просто обратной — к —0(2).
Согласно (11), имеем:
(У* logdet —0(2))* =
= * (—0—1 (2) (—0 (2))) =
= ^ (М0ТиТш—1иМ0М0Т (^ —0(2)) М0) =
= ^ (М'М,^ТиТШ—1иМ'М^ (—0(2))) =
= ^ (Р'иТШ—1иР^—0(2))) =
= tr (5(2) (чк —0(2))).
Аналогично для (У2 logdet —0(2))у.
Таким образом, при г < п + 1 задача (15) в основном сохраняет свой вид. Исключением является лишь неравенство —0(2) ^ 0, заменяющееся на —0(2) ^ 0; соответствующие ему компоненты барьерной функции, ее градиента и обобщенного гессиана находятся согласно (16) и теореме 1.
3.3. РЕПАТРИАЦИЯ РЕШЕНИЙ
Если порядок релаксации к недостаточно велик, исходный метод [5] находит консервативную нижнюю границу значения экстремума. Наиболее разрушительным проявлением аналогичного эффекта в трансформированной задаче является возможность получения финальной конфигурации атомов, ни один из которых не принадлежит допустимой области.8
Чтобы избежать данной проблемы, добавим в трансформированную задачу серию неравенств вида pjС*(ж^) + Л/ ^ 0:
Г
/ * = тш £ pj / (ж ^
* j=l
—0(2) = М^diag(pl,p2,... ,Рг)VТМ0 ^ 0,
Г
(17) —*(2) = £pjGi(жj) ^ 0, г = 1,..., т, j=l
(2) = pjGi(жj) + Л/ ^ 0, г = 1,..., т, ] = 1,..., г,
^2 = £ Рц = 1 ^=1
где Л ^ 0 — параметр, значение которого изначально выбирается достаточно большим для выполнения неравенств (2(0)) ^ 0и
далее систематически уменьшается по мере нахождения промежуточных приближений 2(*) модифицированным методом внут-
8 Например, возможна ситуация, когда в задаче ПН с т > 1 часть атомов удовлетворяет только неравенству ^(ж) > 0, другая часть — только неравенству д2(ж) > 0 и т. д., так что ни один атом не удовлетворяет всем ограничениям одновременно — но при этом неравенства
—*(2) ='П,Г=1 Рзд*(жз-) > 0 выполняются.
9
ренней точки.9
Соответствующие новым неравенствам компоненты барьерной функции logdet —у (2) вычисляются напрямую; компоненты градиента и модифицированного гессиана — с помощью (11):
(18)
(V* -у(2))ж^ = Рц^ (—71 (2) ( 5*(ж,
(у* logdet (2))р^ =^ 1(2)5*(жу ^,
(\72 logdet -у = -р2 ^ (—-1(2) ( аху5*(ж))
1 (2) (5*(жц))) ,
^ 2 (2))ж^кр^ = -рц ^ (—-1(2) (^ (ж))
1(2)5*(жц)) ,
('V2logdet (2))р^р^ = -^ (—- 1(2)5*(ж)—- 1(2)5*(ж));
здесь ж^ и рц в качестве нижних индексов обозначают элементы градиента и гессиана, соответствующие данным переменным; все неуказанные элементы равны 0. Отметим, что неравенства —у (2) ^ 0 добавлены в систему искусственно и не имеют отношения к линейным релаксациям; поэтому использование здесь именно модифицированного гессиана V2 обусловлено не столько намерением построить эквивалентную трансформированную задачу, сколько гарантией положительной полуопределенности соответствующей компоненты Н*.
9 Оптимальная стратегия формирования последовательности зна-
чений Л является предметом отдельного исследования.
4. Примеры
4.1. ЗАДАЧА ПН
Рассмотрим задачу, иллюстрирующую работу нового алгоритма с несвязными областями поиска для т = 1, п = 2:
/ * = шш / (ж) = шт(ж2 + 0,1)2,
Х Х
д1(ж) = 1 - 2ж1 - 2(ж2 - 1)2 ^ 0.
На рис. 1 показана область поиска д1(ж) ^ 0 и отмечены локальные экстремумы / (ж) в данной области, из которых нижний является глобальным.
Решим данную задачу, используя конфигурации с г = 1, г = 2 и г = 3. Выберем случайные начальные позиции атомов в окрестности точки (—0,5; 1). Будем делать серию из 15 шагов с ^ = 1, а после нее — 5 серий по 5 шагов, уменьшая ^ в
4 раза в каждой новой серии. Положим также Л = 1000 (так что компоненты барьерных функций, соответствующие неравенствам
—у (2) ^ 0, не будут оказывать заметного влияния на решение).
Для г = 1, как и в одномерном случае, действие алгоритма эквивалентно локальному поиску с логарифмическими барьерными функциями и ньютоновским выбором направления с использованием модифицированного гессиана (рис. 2).
Для г = 2 и г = 3 отдельные атомы получают возможность переходить из одной компоненты связности в другую за счет временного уменьшения веса, позволяющего матрицам —г(2) оставаться положительно определенными (рис. 3 и рис. 4). Это позволяет одному из атомов найти глобальный минимум; положение остальных атомов впоследствии может стать произвольным, поскольку их веса стремятся к 0.
Дальнейшее увеличение количества атомов в рамках изложенного подхода для данной задачи невозможно, поскольку мы рассматриваем только конфигурации с г ^ п + 1.
Рис. 1. Область поиска
1.5 1.0
0.5 0.0 -0.5 -1.0 -1.5
-1.5-1.0-0.5 0.0 0.5 1.0 1.5 Рис. 2. График х для г = 1
Рис. 3. Графики х и р для г = 2
1.5Г
1.5Ь............................Н
-1.5-1.0-0.5 0.0 0.5 1.0 1.5
Рис. 4. Графики х и р для г = 3
4.2. ЗАДАЧА ПМН
В качестве более реалистичного примера рассмотрим задачу Ж1 из библиотеки COMPleib [10]: нахождение стабилизирующей обратной связи по выходу и = Ку для системы
ж = Аж + Ви,
У = Сж,
где матрицы А, В и С заданы; размерности вектора состояния, управления и выхода равны 4, 2 и 1. Множество решений К = [к1 к2]Т показано на рис. 5.
Модифицируем задачу с целью ее усложнения:
• с помощью замены переменных применим аффинное преобразование к плоскости К для усиления невыпуклости множества решений (что эквивалентно соответствующему изменению матриц А, В и С);
• дополнительно потребуем минимизации к2 + к|; начальные приближения в алгоритме будем генерировать в окрестности допустимой точки, близкой к ложному экстремуму.
Новое множество решений К = [&1 &2]Т показано на рис. 6. Также на графике отмечены локальные экстремумы целевой функции с учетом ограничений. Правый экстремум является глобальным.
В отличие от [5], где данная задача также рассматривалась (в исходной постановке), мы можем напрямую воспользоваться ее представлением в форме ПМН: (А + ВКС)ТР + Р(А + ВКС) < 0, где К е М2х1 и Р е М4х4, Р = РТ > 0, — неизвестные матрицы. Добавив в систему ограничения на неизвестные величины, приходим к следующей задаче ПМН с 12 скалярными неизвестными и 4 неравенствами (левая часть каждого — матрица 4 х 4),
1.0
2.0
0.5
0.0
-10........................<1
-1.0 -0.5 0.0 0.5 1.0
кі
Рис. 6. Решения измененной задачи
0.0
-0.5
-0.5Ь........................А
-1.5-1.0 -0.5 0.0 0.5 1.0
кі
Рис. 5. Решения исходной задачи
из которых одно является билинейным:
Сі(Р, К) = -(А + ВКС)ТР - Р(А + ВКС) > 0, ^(Р,К) = Р - 10-21 > 0,
Сэ(Р,К) = 1021 - Р > 0,
С4(Р, К) = diag(1 + к1, 1 - к1, 1 + к2, 1 - к2) > 0.
Начальные приближения будем формировать из Ко = [-0,9; 0]Т и случайных Ро, находящихся в окрестности решения уравнения (А + ВК0С)ТР + Р(А + ВК0С) = -1. Воспользуемся последовательностью значений аналогичной предыдущему примеру, но через каждые 5 итераций будем делать вспомогательную короткую серию шагов, уменьшая на ней по мере возможности величину Л. Начальное значение последней — Л = 10; в процессе работы алгоритма оно уменьшается до величины порядка 10-3.
Результаты для г = 1, г = 2 и г = 3 приведены на рис. 7, 8
• при г = 1 алгоритм находит ближайший локальный экстре-
• при г = 2 результат, как правило, тот же. Дополнительный атом теряет вес раньше, чем обнаружит правильный экстре-
и 9:
мум;
мум, после чего его траектория становится малопредсказуемой;
• при г = 3 алгоритм достаточно стабильно локализует оба экстремума. Атом в глобальном экстремуме быстро набирает вес, вследствие чего точность нахождения последнего выше, чем точность нахождения ложного минимума.
Таким образом, атомная оптимизация позволяет получать адекватные результаты при решении практических задач с использованием относительно небольших конфигураций атомов.
к\
Рис. 7. График К для г = 1
1.0Г
к\
Рис. 8. Графики К и р для г = 2
к\
Рис. 9. Графики К и р для г = 3
5. Заключение
Мы рассмотрели применение техники атомной оптимизации к невыпуклым задачам на базе полиномиальных матричных неравенств и продемонстрировали следующие ее характеристики.
• Полиномиальная зависимость размера трансформированной задачи от размера исходной задачи ПМ/ПМН и количества атомов позволяет новому алгоритму иметь в целом существенно более низкую вычислительную сложность по сравнению с оригинальным алгоритмом на базе метода моментов за счет отказа от максимально полного исследования области поиска.
• В задачах ПМН интересующего нас вида (охарактеризованного в разделе 3) она позволяет получать адекватный результат при использовании конфигураций атомов относительно небольшого размера. Кроме того, ее одноатомный вариант может использоваться в задачах локальной оптимизации.
• Отметим также, что полученная вычислительная схема применима без изменений в том числе и к задачам, содержащим неполиномиальные функции и неравенства (при
выполнении некоторых базовых требований: в первую очередь, их дифференцируемости и гладкости).
Дальнейшие потенциально перспективные направления исследований включают: поиск оптимального представления и реализации изложенной вычислительной схемы; отработка полученных алгоритмов на существенно нелинейных и невыпуклых задачах теории управления; построение аналогичных вычислительных схем на базе более сложных методов внутренней точки.
Литература
1. ПОЗДЯЕВ В.В. Атомная оптимизация, часть 1: трансформация пространства поиска и одномерные задачи // Управление большими системами. - 2011. - №36. -
С. 39-80.
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. GALEANI S., HENRION D., JACQUEMARD A.,
ZACCARIAN L. Design of Marx generators
as a structured eigenvalue assignment: LAAS-
CNRS Research Report. - 2013. - URL:
http://homepages.laas.fr/henrion/Papers/marx.pdf (дата обращения: 15.05.13).
4. HENRION D., LASSERRE J.-B. Detecting global optimality and extracting solutions in GloptiPoly // Positive polynomials in control. - 2005. - P. 1-18.
5. 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.
6. KORDA M., HENRION D., JONES C.N. Inner
approximations of the region of attraction for polynomial dynamical systems: LAAS-CNRS Research Report. - 2012. -URL: http://homepages.laas.fr/henrion/Papers/roainner.pdf
(дата обращения: 01.12.12).
7. LASSERRE J.-B. Global optimization with polynomials and the problem of moments // SIAM J. on Optimization. - 2001. -Vol. 11, №3. - P. 796-817.
8. LASSERRE J.-B. A new look at nonnegativity on closed sets and polynomial optimization // SIAM J. on Optimization. -2011. - Vol. 21. - P. 864-885.
9. LAURENT M. A comparison of the Sherali-Adams, Lovasz-Schrijver and Lasserre relaxations for 0-1 programming // Mathematics of Operations Research. - 2001. - Vol. 28. -P. 470-496.
10. LEIBFRITZ F. COMPleib: Constraint Matrix-optimization Problem library — a collection of test examples for nonlinear semidefinite programs, control system design and related problems. - 2004. - URL: http://www.compleib.de/ (дата обращения: 01.12.12).
ATOMIC OPTIMIZATION, PART 2: MULTIDIMENSIONAL PROBLEMS AND POLYNOMIAL MATRIX INEQUALITIES
Vladimir Pozdyayev, Arzamas Polytechnical Institute of R. E. Alekseev Nizhny Novgorod State Technical University, Arzamas, Cand.Sc., associate professor (vpozdyayev@gmail.com).
Abstract: We investigate multidimensional optimization problems with polynomial objective function and polynomial matrix inequality constraints and suggest a transformation of the moment-theory-based solution technique. It allows reducing significantly the computational complexity while keeping the ability to solve the problems of the class under consideration.
Keywords: nonlinear programming, matrix inequalities, polynomial inequalities, moment theory.
Статья представлена к публикации членом редакционной коллегии П. С. Щербаковым Поступила в редакцию 09.01.2013. Опубликована 31.05.2013.
