Home    About    Documentation    Install    Newsline    Links    Feedback    Appeal    Guestbook

since 2007/10/24

Отчет и.о. м.н.с. отдела оптимизации ИК НАНУ Крошко Д. Л. о проделанной работе за период июнь-октябрь 2008 г.

1. Реализация блочной схемы Лаптина Ю.П.

В разрабатываемое автором отчета (далее "автором") (на языке Python) ПО OpenOpt была добавлена реализация блочной схемы Лаптина Ю. П. Данная реализация имеет несколько отличий от реализации Лаптина Ю. П., выполненной на С++ с использованием Rational Rose, в первую очередь обусловленных специфическим отличием Python от С++.

Тестирование данной схемы в основном проводилось на нескольких постановках задачи с цепочкой, которая заключается в следующем:
Пусть имеется цепь из звеньев длиной l1, ... lk, связанных шарнирами массой m1, ... mk-1. Если нагрузка на шарнир i превышает Fmax_i, цепь разрушается. Пусть левый конец цепи зафиксирован в точке с координатами (x0, y0). Требуется оптимизировать целевую функцию

1. F = xk -> max (рассматривалась также и постановка xk -> min)
2. F = xk + yk -> max

при ограничениях

a) Только на прочность шарниров: F_i <= Fmax_i
b) На минимальную ординату шарниров: y_i >= Y

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

3. Sum(m_i) -> max (сумма дополнительных нагрузок к шарнирам)

и еще один вариант ограничений:

c) Sum(m_i) = M (либо, как вариант, Sum(m_i) >= M)

Для сравнительного анализа были использованы следующие решатели задач NLP, подключенные автором в настоящий момент к OpenOpt (все приведенные ниже - с открытым кодом и OSI-approved лицензиями):

• IPOPT (Interior Point Optimizer) - наиболее известный бесплатный решатель задач NLP, лицензия CLP, в настоящее время написан на С++, авторы Carl Laird (Carnegie Mellon University) и Andreas Wachter, см. также результаты сравнений IPOPT с коммерческими решателями
• COBYLA (Constrained Optimization By Linear Aproximation) - автор кода (на Fortran) Michael J. D. Powell, лицензия BSD, включен в состав scipy, не использует значения производных
• ALGENCAN (Augmented Lagrangian GENCAN) - лицензия GLP, авторы J. M. Martinez, Ernesto G. Birgin
• slsqp (Sequential Least SQuares Programming) - автор Dieter Kraft, лицензия BSD, включен в состав scipy
• Автор пытался подключить trial-версию коммерческого mosek (у которого недавно появился интерфейс к языку Python) однако, возможно из-за редкой архитектуры (AMD Athlon X2, Linux KUBUNTU) скомпилированная библиотека при вызове выдает системную ошибку обращения к памяти. В будущем автор все же намеревается попробовать ӕто еще раз.
• ralg - реализация р-алгоритма Н.З. Шора, выполненная автором на языке Python с использованием библиотеки numpy
  

Отметим, что некоторыми отличительными особенностями данной реализации являются:

• остановка line search не по признаку тупого угла (что всегда требует вычисления градиента), а по признаку ухудшения точки (как принято в фильтровых методах). Данный подход был одобрен и.о. отдела Стецюком П. И., который высказал намерение использовать его (в несколько ином контексте и комбинации с другой идеей) в фортран-версии р-алгоритма.
• отсутствие вычисления градиента целевой функции в infeasible-точках (там, где есть активные ограничения). Расчет значения целевой функции в них производится только для вывода текстовой и/или графической информации между итерациями, а также чтобы дать возможность проверить некоторые критерии остановки, определенные пользователем и/или ядром OpenOpt. В отличии от схемы, описанной в статье "Использование модификации р-алгоритма для нахождения глобального минимума полиномиальных функций", Шор Н. З., Стецюк П. И. "Кибернетика и системный анализ" № 4/1997, стр 36, формула (16), для вычисления направления растяжения используется следующая вектор:
  • разность градиентов целевой функции, если обе точки текущей и прошлой итерации не имеют активных ограничений
  • разность градиентов максимальных ограничений в точках текущей и прошлой итераций, если обе точки имеют активное ограничение
  • градиент максимального активного ограничения в одной из ӕтих точек (текущей и прошлой итерации), если одна и только одна точка не имеет активных ограничений.

Следует обратить внимание, что

• Данный подход является ӕквивалентом целевой функции F(x) = f(x) + S*max(c_i(x)) при очень большом коӕффициенте штрафа S; в данных обозначениях f - целевая функция, х - вектор оптимизируемых переменных, c_i - ограничения вида c_i(x) <= 0.
• В настоящее время ограничения на равенства h_i(x) = 0 обрабатываются как -contol <= h_i(x) <= contol, где contol - требуемая точность (по умолчанию в OpenOpt используется 10^-6).
• Данная реализация снижает число умножений р-алгоритма на каждой итерации с 5*n² до 2*n² + 3*s*n, где s - число ненулевых ӕлементов градиента максимального ограничения (для feasible-точек - целевой функции).
  
• Несмотря на ряд тестов, в которых данная реализация р-алгоритма превосходит остальные протестированные решатели, она все еще имеет ряд недостатков, работа над устранением которых продолжается. Одним из них является неприспособленность классической для р-алгоритма регулировки шага к переходам между feasible и infeasible итерационным точкам (в openopt-реализации), т.к. модуль градиента целевой ф-ции и максимального ограничения (в соседних итерационных точках) могут значительно отличаться, что приводит к плохому выбору стартового шага для поиска по направлению.
• Все протестированные решатели (кроме ralg) требуют, чтобы целевая функция в стартовой точке была определена.
  

Перед тем как привести ниже результаты численных ӕксперементов, сделаем следующую краткую характеристику использованных решателей исходя из результатов их тестирования на блочных задачах без доопределения ОДЗ нелинейных функций:

• COBYLA - ӕффективен только для задач небольшой размерности, на бОльших чем nVars ~ 10 сильно уступает по скорости
• slsqp - в ряде случаев показывает очень большое быстродействие (быстрее всех остальных решателей в 2-3 раза), однако работает нестабильно. Часто останавливается с кодом "LSQ matrix is singular to working precession" (впрочем, достаточно близко к оптимальной точке), что, возможно, обусловлено ошибками разработчиков scipy или f2py его при стыковке с Python. Поӕтому в ряде случаев может служить хорошим средством для поиска приблизительного решения, чтобы потом попытаться дорешать его другим решателем. Однако на задачах, где целевая функция и/или нелинейные ограничения могут возвращать NaN (not a number), ӕффективность slsqp оказалась очень низкой, а процент решенных задач очень небольшой.
• ALGENCAN - несмотря на декларирование о хорошем механизме обработки областей определения нелинейных функций (что отображено даже в его документации и интерфейсе посредством флагов), показал довольно низкую ӕффективность на протестированных блочных задачах, как впрочем и на ряде обычных, где функции заданы на всем R^nVars. Ӕто относится как к скорости, так и к проценту решенных задач. Следует отметить, что в ряде случаев ALGENCAN 1.0 показывал лучшие результаты. Разработчики ALGENCAN продолжают работу над текущей версией 2.0.5beta и, возможно, ӕта ситуация будет улучшена.
• IPOPT - де-факто самый известный из бесплатных (и, более того, с OSI-одобренной лицензией) решатель задач NLP на протестированных задачах, наряду с ralg, показывал наиболее стабильную и быструю работу. Была использована версия 3.3.5, в то время как недавно уже вышла более новая - 3.5.4. Однако в ряде случаев и у него наблюдались либо остановки достаточно далеко от оптимальной точки, либо очень пологие плато на траектории сходимости, за которыми иногда следовал скачек к оптимальной точке, что вообще говоря свидетельствует о ненадежности, т.к. критерии остановки могут сработать до ӕтого скачка. Из 3-х решателей систем линейных уравнений, которые IPOPT требует при установке, был выбран MUMPS, т.к. только он распространяется с OSI-одобренной лицензией (GPL).
• ralg - по стабильности показал примерно те же результаты, что и IPOPT; по скорости отличие в среднем состояло в 2-3 раза, иногда как в одну сторону, так и в другую (следует также иметь в виду, что Python, даже с учетом использования numpy, все же уступает по скорости С++). Ситуация с плато на траектории у ralg практически никогда не возникала (на протестированных блочных задачах с цепочкой).

Сдедует отметить, что все способы доопределения, предлагаемые Лаптиным Ю. П., приводят к негладкой задаче, соответственно от гладких решателей не следует ожидать хорошего и стабильного решения. Однако наряду со специализированным ralg, достаточно хорошие результаты показал IPOPT, а также (для небольших nVars) - COBYLA. Вместе с тем, автору не удалось получить ни одного примера, где доопределение, которое проводилось по правилу F(x) = f(x) + S * distance(x, dom f), где S - величина штрафа (были опробованы несколько достаточно больших значений), ускорило или улучшило работу какого-либо из решателей.

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

(to be continued) - я в настоящее время, перед тем как печатать окончательные результаты сравнения, пытаюсь в ралге подкрутить регулировку шага для старта в поиске по направлению.

2. Дальнейшие планы по развитию блочной схемы (на английском доступны по ӕтому адресу):

• автоматическая проверка производных для блоков (в настоящее время проверка доступна только для классического синтаксиса задач, см пример)
• удобная обработка фиксированных переменных
• проверка идеи о растяжении пространства внутри каждого блока

3. Продвижения в разработке OpenOpt
Наряду с работой над реализацией блочной схемы Лаптина Ю. П. за отчетный период автором была проделана следующая работа над дальнейшим улучшением OpenOpt:

• Некоторые изменения к ralg (особенно связанные с обработкой линейных ограничений Ax <= b, Aeq x = beq, lb <= x <= ub)
• Подключен ALGENCAN 2.0.3 (вместо ALGENCAN 1.0)
• Исправлен вывод графики для задач NLSP (системы нелинейных равенств) с ограничениями
• Добавлен параметр шкалирования в решатель lpSolve
• Создан новый класс LLAVP (linear least absolute values)
• Улучшена обработка нелинейных функций с ограниченной областью определения
• Решатель galileo теперь может обрабатывать задачи с целочисленными переменными
• Подключен еще один решатель для класса GLP (поиск глобального ӕкстремума): pswarm
• Добавлены следующие преобразователи (задач из одного класса в другой): llsp2nlp, llavp2nsp, llsp2qp, llsp2nlp, lsp2nlp, qp2nlp, lp2nlp, nlsp2nlp
• Добавлена удобная возможность обработки задач максимизации через p.goal = 'max' или 'maximum', соответствующие изменения в коде проведены для текстового и графического вывода.
• Несколько исправлений и чистка кода.

(to be continued)