Единый параллельный алгоритм методов приведенных направлений для решения задачи нелинейной оптимизации с ограничениями-неравенствами Текст научной статьи по специальности «Математика»

УДК 004.4.032.24:519.853

О. В. Жеребцова

ЕДИНЫЙ ПАРАЛЛЕЛЬНЫЙ АЛГОРИТМ МЕТОДОВ ПРИВЕДЕННЫХ НАПРАВЛЕНИЙ ДЛЯ РЕШЕНИЯ ЗАДАЧИ НЕЛИНЕЙНОЙ ОПТИМИЗАЦИИ С ОГРАНИЧЕНИЯМИ-НЕРАВЕНСТВАМИ

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

Введение

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

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

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

1 Методы приведенных направлений для решения нелинейных задач с ограничениями в виде неравенств

Приведем краткое описание теоретически обоснованных в [2-6] методов приведенных направлений для решения задач оптимизации с ограничениями-неравенствами.

Известия высших учебных заведений. Поволжский регион Рассматривается задача

шш / (х); 0 = | х eEn | (x )<0, у еJ}, (1)

здесь Еп - п -мерное евклидово пространство; J =|1, т} - конечное множество индексов; функции fj (х), у е J и|0} достаточно гладкие.

Под решением задачи (1) понимается точка х* , удовлетворяющая необходимым условиям минимума первого порядка:

Л (х*)+ £ А*// (х ) = 0

jеJ

ЗА* > 0: ) (х* ) = 0, (2)

Л) (х* )< 0, уе J.

При этом предполагается, что в каждой точке х* , являющейся решением задачи (1), выполняется условие строгой дополняющей нежесткости, т.е. в (2):

А* > 0, уе J*, J* =|уе J|fj. (х*) = 0}.

Для решения задачи (1) строится итерационный процесс

хи+1 = хк + s(^), к = 0, 1, 2,...; ( = 1, 2,

где s(()(^к) - траектория движения из текущей итерационной точки хи; tк > 0 - шаг вдоль траектории.

В соответствии с [1] траектории движения к новой итерационной точке определяются как

51 () = ts - линейная;

52 () = ts - 0,5^Яц - криволинейная,

где 4 = (*))^еJ(х,в).

Для построения приведенного направления 5 в текущей итерационной точке используются варианты построения множества активных ограничений J(х,е), е> 0, способы выбора которых для соответствующих методов приведены в [2-7], при этом предполагается, что градиенты (х), у е J(х,е)

являются линейно независимыми для всех х е Еп .

Направление движения к следующей итерационной точке определяется следующим образом:

5 = -Р(х)г-Я(х)( + /(х)). (3)

Векторы г є Еп г , V є Ег , г = | 3 (х, е) |< п - параметры направления, /(х) = [/] (х))^(х е). Матрщы Р(х) и Я(х), размерности пх(п - г) и п х г соответственно определяются условиями:

ЛТр(х) = О, АТЕ(х) = 1г, А = ((х)) ,є_,(хе) (4)

В работах [1-7] предлагается построение матриц Р (х), Я (х), удовле-

Т

творяющих условиям (4), основанное на LQ -разложении матрицы А :

Р(х) = И$, Я(х) = Н^-1,

АТН = (^0) 0, 0Т = (, 02). (5)

Т

где О = НН - симметричная положительно определенная п х п матрица; Н -

нижняя треугольная, L - левая треугольная г х г матрицы; 0 - ортогональ-

ная п х п матрица.

Согласно (3), построение приведенного направления 5 в текущей итерационной точке сводится к построению матриц Р, Я и векторов г є Еп г , V є Ег .

Следующие формы векторов г, V обеспечивают убывание функции

/о (х) при движении вдоль траектории 5(і) (ґ), і = 1,2 :

г = -£ (х); (6)

г = -(Р(х)Т° Р(х)) Р(х)Т/0(х),

О =

/о (х)+ 2 и]/] (х)

Іє3 (х,е)

(7)

V = [-и (х)- /(х)]+ , (8)

у хх

+

здесь

8(х) = -Р(х) /0(х), и(х) = -Я(х) /0(х). (9)

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

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

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

а) способа построения множества активных ограничений;

б) расчета функции выигрыша и штрафного параметра;

в) способа вычисления параметров, определяющих направление,

что позволяет осуществить единообразную программную реализацию описанных методов.

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

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

2 Параллельный вычислительный алгоритм методов приведенных направлений

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

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

Рис. 1 Блок-схема единого вычислительного алгоритма методов приведенных направлений

0. НАЧАЛЬНАЯ ПРОГОНКА. Вычисления значений функций задачи в начальной точке (рис. 1, шаг 0.2) независимы и могут производиться параллельно. Для нахождения множества индексов активных ограничений (рис. 1, шаг 0.3) все сравнения тоже могут выполняться одновременно.

1. ПОСТРОЕНИЕ НАПРАВЛЕНИЯ ДВИЖЕНИЯ. На первом шаге

(рис 1, шаг 1.1) сформируем матрицу А = (/,-{х)) состоящую из

' 'jeJ(х, е)

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

Т

матрицу А функций по переменной х/ . Тогда матрица А распределена по столбцам, их количество не зависит от размера множества активных ограничений J {х, е) и равно п:

АТ =

' Л >

Эх1

Эх2

дхп

V п /

В случае реализации ЬО -разложения матрицы А (рис. 1, шаг 1.2) для последовательных компьютеров выбор алгоритма выполнения разложения не существенен. Однако для параллельной реализации выбор того или иного алгоритма выполнения ЬЦ -разложения обусловлен: а) способом хранения исходных данных на момент начала работы алгоритма; б) способом хранения результирующих данных. Кроме того, т.к. ЬЦ -разложение является частью

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

Т

Процесс ЬЦ разложения матрицы А можно представить как после-

ТТ довательное умножение матрицы А справа на ортогональные матрицы О ,

которые являются преобразованиями Хаусхолдера:

Ь = АТЦТ = АТОС02 •...• еТ_ь 2 = 6122 •...• 0,-1;

Ь{0)= АТ; Ь{к+1 = Ь^О+ь к = 0, 1, 2, ..., п _ 2.

Преобразованием Хаусхолдера называется матрица вида

Т Т I _ ww , ^ ^ = 2,

где w - вещественный вектор-столбец. Матрица Хаусхолдера симметрична и ортогональна [9].

Полагая

Т wk

1к =10, •••, 0, 1кк С, 1к,к+1, •••, 1к,п

(10)

где

С = _

Е ( м2

j=к

ц =

С2 _ С I

(к) к к

(11)

получаем

(к) к ,к

+1 ((к) )2

І=к+1

= ц

х ( кк)

І=к

~2Ы й +с2

л

= 2Ис2-с/ Ц ) = 2.

Т Т

Тогда Ок+1 = I _ WkWk является матрицей Хаусхолдера. Применяя пре-

образование Хаусхолдера на к -м шаге, имеем

Ь(к+1) = Ь(к) е[+1 = Ь(к) _ Ь(к) WkWС .

Учитывая, что

\2

(12)

‘к,к + X (

І=к+г

= ^ ( " $) = ^ ’

(к)

7 (к)

где ¡к - соответствующая строка матрицы Ь ', имеем

= с.

¡к“’=¡{к _^( _ с )4" ч

¡к‘+1)= -Л* ¡{] ¡кк)“> = 0, ./ = к +1, .... п.

Таким образом, в результате применения преобразования Хаусхолдера на к -м шаге получается новая матрица с нулями в позициях (к, /), где j = к + 1, ..., п, при этом нули, полученные на предыдущих шагах, сохраняются в соответствующих позициях.

Т

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

Рассмотрим к -й шаг алгоритма ЬО -разложения методом Хаусхолдера:

I. Найдем значения с и * из (11). Для этого параллельно вычисляются

значения {{у ) , У = к, . .,п . Находится сумма Е (к/ ) операцией сдваи-

/=к

вания и результат рассылается всем участникам вычислений. Рассылается значение элемента ¡кк и вычисляются значения с и * . С учетом (10) параллельно формируются компоненты вектора Wk .

II. Выполним преобразования /-й строки. Вычисляются параллельно

п

dj ¡/к, / = 1,2,..., п, и операцией сдваивания находится и = Е^, ре-

/=1

зультат рассылается всем участникам процесса вычислений. Преобразуем /-ю

,(к+1) ,(к) Т

строку: ¿> ' = ¡1 ' - и w .

к

III. Аналогичным образом получаем текущую матрицу 2 = 61 •... • Ок.

IV. Повторяем шаги П-Ш для / = к +1,..., п .

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

а) на шаге I потребуется 3р пересылок;

б) на шаге II - 2р пересылок.

Итого, с учетом повторений цикла, имеем 3пр + п(п + 1)р пересылок. Проведем грубую оценку эффективности алгоритма Ь2 -разложения по методу Хаусхолдера. При последовательном выполнении вычислений прово-

3 2

дится примерно 2п + 5п операций (слагаемые, содержащие степени п, меньшие 2, опустим). Тогда, считая, что ^ - время выполнения одной операции, (у> 1) - время, затрачиваемое на пересылку одного элемента векто-

ра или матрицы, получим следующую оценку эффективности:

2п3+ 5п2

Ер =-------------------;------Г .

2п3 + 5п2 + у р2 (2 + 4п)

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

По окончании ЬО -разложения результирующие матрицы Я и Р размерностей п х г и п X (п - г) соответственно хранятся по столбцам. Будем иметь матрицу (Я;Р) размера пхп .

Для нахождения векторов и(х) = -ЯТ/¿(х) и g(х) = -РТ/¿(х) (рис. 1,

шаг 1.3) выполним матрично-векторное умножение -(Я;Р) /)(х) методом скалярных произведений [9]. Тогда первые г независимых скалярных произведений ^-г/,/о(х)^, / = 1,..., г сформируют компоненты вектора и (х), а

оставшиеся (п - г) произведений ^_p¡,/)(х)^, ¡ = 1,..., п - г - компоненты вектора g (х). Здесь г/, р1 - соответствующие столбцы матрицы (Я;Р). Итоговый вектор (и(х); g(х) ) будет распределен между процессами покоординатно.

Для нахождения направления движения 5 из (7) формируется вектор

Вычисление произведений под знаком суммы является независимой операцией и может выполняться одновременно. Результатом каждого умножения является вектор-составляющая направления, которые после суммирования методом сдваивания дают итоговый вектор 5 .

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

2. ВЫЧИСЛЕНИЕ ДЛИНЫ ШАГА. Для нахождения длины шага (рис. 1, шаг 2.1) вычислим новую точку с текущим значением шага, вычисления дублируются необходимое число раз. Следующие шаги 2.2-2.4 аналогичны шагам 0.2-0.4 начальной прогонки, в дальнейшем полученные значения будут использоваться на этапе нахождения направления движения в следующей точке.

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

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

1. Ижуткин, В. С. Методы приведенных направлений для задачи нелинейного программирования / В. С. Ижуткин, М. Ю. Кокурин // ЖВМиМФ. - 1988. - Т. 28. -№ 2. - С. 1799-1814.

2. Ижуткин, В. С. Методы приведенных направлений с допустимыми точками для задачи нелинейного программирования / В. С. Ижуткин, М. Ю. Кокурин // ЖВМиМФ. - 1990. - Т. 30. - № 2. - С. 217-230.

3. Ижуткин, В. С. Методы приведенных направлений на основе дифференцируемой штрафной функции для задачи нелинейного программирования /

( + (1 ((,..., (г + /г (х)) , 21,..., 2п-г,) длины п, и методом линейных комбинаций [9] выполняется умножение (рис. 1, шаг 1.4):

J г=1 /=1

I=1

Заключение

Список литературы

B. С. Ижуткин, М. В. Петропавловский // Известия вузов. Математика. - 1994. -№ 12. - С. 50-59.

4. Ижуткин, В. С. Методы приведенных направлений на основе модифицированной функции Лагранжа для задачи нелинейного программирования / В. С. Ижуткин, М. В. Петропавловский // Известия вузов. Математика. - 1995. - № 12. - С. 33-42.

5. Ижуткин, В. С. Методы барьерных функций и центров для задачи нелинейного программирования / А. В. Блинов, В. С. Ижуткин, М. В. Петропавловский // Известия вузов. Математика. - 1996. - № 12. - С. 32-41.

6. Ижуткин, В. С. Метод внешних центров с использованием приведенного градиента для задачи нелинейного программирования / А. В. Блинов, В. С. Ижуткин // Известия вузов. Математика. - 1998. - № 12. - С. 49-57.

7. Ижуткин, В. С. Комбинированный метод приведенного градиента для задачи нелинейного программирования и его программная реализация / О. В. Бастракова (О. В. Жеребцова), В. С. Ижуткин // Математическое программирование и приложения : материалы XI Всероссийской конференции. - Екатеринбург, 1999. -

C. 131-132.

8. Ижуткин, В. С. Мультистадийный метод приведенных направлений для задачи нелинейного программирования / О. В. Бастракова (О. В. Жеребцова), В. С. Ижут-кин // Распределенные системы: оптимизация и приложения в экономике и науках

об окружающей среде : сборник докладов Международной конференции. - Екатеринбург, 2000. - С. 139-141.

9. Ор тега, Дж. Введение в параллельные и векторные методы решения линейных систем / Дж. Ортега. - М. : Мир, 1991. - 366 с.