Опыт декомпозиции метода конечных элементов с использованием теории структурированных оптимизационных задач Текст научной статьи по специальности «Математика»

Опыт декомпозиции метода конечных элементов с использованием теории структурированных оптимизационных задач

А. С. Величко (vandre@dvo.ru), Е. А. Нурминский (nurmi@dvo.ru) Институт Автоматики и Процессов Управления ДВО РАН

Аннотация

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

Введение

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

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

Оценивая состояние теории и практики методов декомпозиции можно отметить, что особенно полно развита теория декомпозиции задач линейного программирования, где ее применение может существенно уменьшить вычислительные затраты [6], [7]. В алгоритмах типа симплекс-метода для решения задач линейного программирования с п переменными и т ограничениями пессемистические оценки обьема вычислений имеют экспоненциальный порядок вида 0(ап+т) стандартных операций, где а > 1, В последние годы большое внимание уделялось развитию несимплексных алгоритмов с полиномиальной сложностью, для которых такие оценки имеют вид 0(пк), где показатель степени обычно порядка 5, Однако, эта оценка предсказывает для практических задач5 где тъ ^ 10 , также весьма значительный объем вычислений. Кроме этого в этих алгоритмах существенным являются и требования к памяти ЭВМ (порядка иг).

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

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

1 Постановка задачи

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

где г а, -/; векторы переменных задачи размерностей Ма и Мв соответственно, .г -вектор связывающих переменных задачи размерности са, <'н векторы длины Ма и Мв соответственно, А а, Ав, В а, Вв - матрицы соответствующих размерностей, <1 л -вектор длины К а, <1ц вектор длины Кв.

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

ппп сАгА + свгв а а* а + вах < йа

Авгв Ввх ^ <1 ц -л > 0. > 0. х > 0,

(1) (2)

(3)

(4)

Без выделения блоков задача ( 1 )-( 1 ) запишется как minez при Az < <1. z > 0, где z = (zA,zB,x), с=(сА,св, 0), (I = (>/л-'//;)•

4 = f ° Вл \ 0 Ав вв

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

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

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

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

Введя функции

fA(x) = niinr л: л. fB(x) = min cBzB, (5)

ZA '¿В

AaZa < da - bAx Aazb < (ii! - bBx

> 0, x > 0 zB > 0, x > 0

получим эквивалентную ( 1 )-( 1) задачу:

mm{fA{x) + fB{x)}. (6)

Как известно, fA(х) и fB(x) - выпуклые и кусочно-линейные функции [4]. Введем для них сопряженные функции hA(p) и hB(p):

Ь>а(р) = Я(ж) = maх{рх - fA(x)}, hB{p) = f *B{x) = тах{рх - fB(x)}.

Отсюда в частности следует, что

~Ьа{~р) = rmn{cAzA+px}, - hB{p) = min{cBzB - рх}, Aaza + Вах < dA ABzB + Ввх < dB

> 0, х > 0 zB > 0, х > 0

то есть вычисление hA(^p), hB(p) также сводится к решению задач линейного программирования, но уже по всем переменным с дополнительными стоимостями для связывающих переменных в качестве параметров.

Задачу (6) можно легко переписать в терминах сопряженных функций:

min(fA(x) + fB(x)} = min {¡А(хг) + fB(x2)} =

x x1—x2=Q

= min max{fA(xl) + pxl — px2 + fB(x2)} =

х1—х2=0 P

= max{min(//4(a;) + px) + min (fB(x) — px)} = = max{~ max(^pa; — fA(xj) — max(px — fB(x))} = = min{max(-pa; - fA(x)) + max(px - fB(p))} = = - min{hA(-p) + hB(p)}.

Последнюю задачу в этой цепочке равенств

min{hA(-p) + hB(p)} (7)

будем называть двойственной к (6),

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

2 Алгоритм двойственных отсечений (АДО)

Следуя [4] сформулируем алгоритм двойственных отсечений (АДО) для задачи (6), Суть этого подхода заключается в замене некоторых функций в (6)-(7) на их аппроксимации, получаемые в процессе вычислений и обмена, В (6) аппроксимируется функция /В(х), в (7) - сопряженная функция другого блока кА{р). Эти аппроксимации будем обозначать как ¡в{х) и ккА(р) соответственно. Собственно алгоритм формулируется следующим образом:

Шаг 1 Выбрать начальные приближения х°, р°.

Шаг 2 Для начального приближения р° вычислить Ь,в(р°). Определить начальные аппроксимации ]'п 1 (>) = —оо, //л' (}>) = —оо. Положить к = 0,

Шаг 3 ЕСЛИ текущая точка хк не оптимальна, то выполнять: BEGIN

Решить задачу

1в(хк) = min cBzB (8)

ZB,X

ABzB -Ь Ввх fi dB h

х — X

zB > 0, ./• > 0

и определить рк е д/в(хк), ¡гв(рк), Модифицировать прямую аппроксимацию 1В(Х) = тах{/£_1(а;), ркх - Нв(рк), ркх - Нв(рк)}, Решить задачу

ппп{Ых) + ¡к(х)} (9)

и определить хк+1, /А(хк+1), Решить задачу

^ЬА(^рк)= тт(сАгА + ркх) (10)

АА2А + ВАх < йА zA>0, х>0

и определить хк е дкА(^рк), /А(хк), Модифицировать двойственную аппроксимацию ¡ъкА(р) = тах{/1д_1(р), рхк — /А(хк), рхк+1 — /А(хк+1)}, Решить задачу

mm{hkA(-p) + hB(p)} (И)

и определить pk+1, hB(pk+1), • к = к + 1 • Перейти к шагу 3.

END

ИНАЧЕ хк - решение задачи (6).

3 Вычислительные аспекты АДО

Важными вычислительными задачами в АДО являются вычисления субградиентов функций /в(-); Ь'а(') в (8) и (10). Из стандартной теории двойственности линейного программирования [5] легко получить следующее утверждение, которое мы все же докажем в целях самодостаточности работы.

Пусть fB(x) определена соотношением (5). Тогда для некоторого фиксированного вектора х ее можно представить в эквивалентном виде:

fß{x) = min cBzB (12)

ZB ,x

ABzB -Ь Ввх fi (Iи

> 0. X > 0.

Теорема 3.1 Пусть задача (12) разрешима для заданного х, и р* является вектором оптимальных двойственных переменных, соответствующих ограничению х = х в задаче (12). Тогда р* G dfB(x), и обратно, любой вектор р* G dfB(x) является вектором, оптимальных двойственных переменных, соответствующих ограничению х = х в задаче (12).

Доказательство:

Покажем сначала, что если р*в является вектором оптимальных двойственных переменных, соответствующих ограничению х = х в задаче (12), то р*в G dfB(x). В силу теорем двойственности для х G dorn fB

fB(x) = maxdu + xpB (13)

«> Рв

< с

и'Вв + рв < 0 и < 0.

Тогда

f(x) — f(x) > du* + хр*в — (du* + хр*в) = (x — x)p*B,

и, следовательно, p*B G dfB(x).

Докажем теперь обратное утверждение. Пусть р*в G dfB(x). Известно, что р*в G д/в(х) тогда и только тогда, когда

1в(х)^хр*в = -кв{рв),

где hB - функция, сопряженная для fB. По определению

-^в (рв) = min cBzB - р*вх

ZB уХ

AßZß ~Ь BßX dß ZB> o, > o.

для которой двойственная задача имеет вид

Ь>в(Рв) = та хс1и

и

и'Ав < с и'Вв < ^р*в

(14)

и < 0.

Сравнивая задачи (13) и (14), легко заметить, что при фиксированном р*в выполняется —Нв{р*в) = 1в{х) — хр*в, ограничения этих задач совпадают, следовательно, р*в является вектором оптимальных переменных в задаче (13) и, следовательно, вектором двойственных переменных, соответствующих ограничению х = х в задаче (12). ■ Вычисление Нв(рк) в алгоритме двойственных отсечений приводит к задаче

Из (9) ясно, что если (гд, хк, V*) - оптимальное решение задачи, то /д(хк) = сдгкА. В силу теоремы 3.1 нахождение хк е д}гд(—рк) приводит к вычислению оптимального вектора хк в задаче

Ьв{рк) = - тт(свгв - ркх) = ркхк - свгкв

Авгв ~Ь Ввх ^ (1ц > 0. .г > 0.

где хк) - решение задачи (12).

Задача (9) может быть представлена в следующем виде:

ЛАгА + ВдХ < йд

V > ргх — ¡1В{р1), 0 < i < к

V > ргх — ¡1В{р1), 0 < i < к

ZA> 0, х > 0.

Ьл{-Рк)= тт{сАгд+ркх)

(15)

Ал^л + ВдХ < ¿А гА > 0, х > 0.

По определению функции /а(х)

¡д(хк)= ПППГ л; 1

Адгд + ВАхк < с!д гА>0

Аа^а < <1л - ВАх-

гА>0-

пнисАгА

.к

Решение задачи (11) не позволяет сразу получать переменные Zß иг, поэтому запишем эту задачу в терминах сопряженной переменной х.

~ + hB(p)} = ~ min Ш{р1) + hB(p2)} =

Р pL+p2=0

= — min таx{hkÄ(pl) — plx — p2x + hB(p2)} =

p1+p2=Q x

= — max{min(/i^(p) — px) + min (hB(p) — px)} = = — max{~ таx(px — hkÄ(p)) — таx(px — hB(p))} = = min{fB(x) + таx(px - hkA{p))}.

Последнюю задачу на максимум представим как

maxim; — v)

p,v

ev > p'Xk - Fk ev > p'Xk - Fk,

где Xk = {x°,...,xk), Xk = (x°,...,xk), Fk = {f{x°),...,f{xk)), Fk = (f(x°),...J(xk)).

Перейдя от переменных (p, v) к вектору двойственных переменных Л, получим двойственную задачу вида

к-

min FaA

а л

•г - .V A = 0

еА = 1, Л > 0,

где Хк = (Хк,Хк), Fk = (Fk,Fk).

Окончательно получаем эквивалентную для (11) задачу:

min(cbzb + F%\) (16)

ZB, X, а

AbZ в И н-1 < dB (17)

x - XkX = 0 (18)

еА = 1 (19)

zB> 0, а; > 0, А > 0. (20)

Вектор pk определяется как оптимальные двойственные переменные ко второй группе ограничений задачи. Ясно, что если (z|, хк, А*) - оптимальное решение задачи, то hß(pk) = Ркхк — cBZß, поскольку функционалы сопряженных задач равны с точностью до знака.

Заметим, что вычисление рк G dfB(xk) представляет собой частный случай задачи (16)-(20), когда переменная Aj, соответствующая хк, равна 1.

В приведенных выше формулировках /д, fB являются выпуклыми кусочно-линей-11ы.\ш функциями, и при этом можно доказать конечную сходимость алгоритма [4].

Теорема 3.2 Пусть fA, fB являются конечными, выпуклыми и кусочно-линейными, все подзадачи АДО являются разрешимыми для любой итерации к. Тогда задача (1)-(4) имеет конечное решение, и для некоторого к

min{fA(x) + fB(x)} = fA(xk) + fB(xk).

Доказательство:

По построению последовательность {/¿(ж)} монотонна и ограничена сверху:

Ш < fkb+1(x) < fB{x) < +оо.

С другой стороны

fB(xk) < max{fk(xk), fB(xk)} = fB+l(xk). Отсюда fB(xk) = fB+l(xk). Тогда

h(xk) + fB(xk) = fA(xk) + fk+l(xk) = 1шп{/л(а;) + fkB+l(x)} < min{fA(x) + fB(x)},

откуда следует, что xk решает задачу min{/^(a;) + fB(x)}.U

В алгоритме мы не конкретизировали условие прекращения его работы. В качестве признака оптимальности полученного решения можно проверять отделяемость приближенного решения (z¿, zkB.; хк) на к-ом шаге алгоритма вновь добавляемым отсечением.

4 Модификация АДО для задач с нетривиальным носителем

Во многих постановках область изменения связывающих переменных неявно ограничена требованиями допустимости задач вычисления функций fA, hB. Особенно явно это проявляется в декомпозиции задач решения систем уравнений, где подзадачи неминуемо оказываются либо пере- либо недоопределенными. Это приводит к неявно заданным нетривиальным ограничениям на связывающие переменные, которые в предыдущей версии алгоритма [4] не учитывались. Этот эффект весьма существенен, на практике во время работы алгоритма, особенно на начальных итерациях, вероятно, а фактически неизбежно, возникновение недопустимых подзадач (8), (10). Характер неразрешимости этих задач может быть двоякой природы: допустимое множество задачи может оказаться пустым, либо ее целевая функция не ограничена на допустимом множестве задачи, то есть допустимое множество задачи содержит направление бесконечного убывания целевой функции (для задачи минимизации).

В АДО вычисление рк G dfB(xk) может привести к задаче с пустым допустимым множеством, а вычисление хк G Oh }>''') - к задаче, где целевая функция не является ограниченной. Ясно, что в этих случаях становятся неопределенными величины hB(pk), 1л(хк). С точки зрения теории двойственности случай неограниченной задачи может быть сведен к вопросу о коррекции алгоритма в случае возникновения подзадачи с пустым допустимым множеством на некотором шаге, поэтому рассмотрим сначала первый случай.

Решить проблему неразрешимости подзадачи алгоритма можно добавляя дополнительные отсечения, которые отделяют область недопустимых аргументов хк в задаче рк G OfB(xk) от допустимых векторов, множество которых описывается условием Íb(x) < +00, //;(•'•) / ^00, Такие отсечения отделяют область недопустимых аргументов хк от эффективного множества функции fB(x), которое обозначим как dorn fB.

Определение 4.1 Гиперплоскость рх = а назовем отделяющей точку х от множества Ф; если рх < а < рх для х G Ф. Такую гиперплоскость будем также называть отсечением,.

Множество dorn / описывается с помощью задачи поиска допустимого базиса в задаче линейного программирования.

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

Положим ||sj^ = X) |s¿| и запишем задачу поиска допустимого базиса при нахождении

pk G dfB(xk) в виде

ф(хк) = min РЦ, > 0

(21) (22)

(23)

(24)

ABzB Ввх ^ <1 ц

X + s = хк

> 0. X >0.

Функция ф(х) конечна и выпукла и ф(х) = 0 для х G dorn fB. При qk G дф(хк)

qk(x - хк) < ф(х) - ф(хк).

Если х G dorn fB, то

ф{хк).

(25)

Поскольку ф(хк) > О, то

qkx < qkxk - ф(хк) < qkxk

для х G dorn fB. Отсюда следует, что (25) - отсечение, отделяющее хк от dorn fB. Задача отделения точки хк от dorn fB, как правило, имеет множество решений. По теореме 3,1 qk находится как двойственный вектор к ограничению (23) в задаче поиска допустимого базиса. Однако, qk, построенный как двойственный вектор к ограничению (23) в задаче (21)-(24), обладает специальным свойством оптимальности.

Определение 4.2 Назовем гиперплоскость в определении 4-1 плотным отсечением, точки х от множ.ества dorn/, если дополнительно выполнено

а = min «i = max рх " I X G dorn /

рх < «i X G dorn /

Рассмотрим отсечение (25), qk G дф(хк), и

а* = max qkx X G dorn fB

Рис, 1: Отсечения, отделяющие вектор хк от dom/в, Пунктирная линия - плотное отсечение.

Теорема 4.1 Отсечение (25) является плотным, т.е. а* = qkxk — ф(хк).

Д о к а з а т о л ь с т в о :

Двойственной задачей для (21)-(24) будет задача линейного программирования вида

ф(хк) = та х(с1,и + рхк)

и,р

и'Лв < О ь'Вв+р < О

^е < р < е V < О,

где е - вектор из единиц соответствующей размерности.

Обозначая ее решение через (</''. й). имеем ф(хк) = дкхк + с1й. Так как дк е дф(хк)

■■ —г1й или

и

(х) — ф(хк) > qk(x — хк), то для х G dorn fB имеем qkx < qkxk

яг

qkx + du < О,

Для того, чтобы показать плотность данного отсечения, достаточно показать, что

max (qkx + du) = О, х G dorn fB

Но если х G dorn fB, то ф(х) = 0, тогда max (qkx + dü) = 0, ■

х G dorn fß

Отсечение (25), аппроксимирующее dom/ß, добавляется к ограничениям задачи (9), которые препятствуют возникновению недопустимости на следующих итерациях АДО, После накопления отсечений допустимости областью изменения связанных переменных будет множество, показанное на рис, 2,

Рис, 2: Плотные отсечения, аппроксимирующие dom fB.

Вычисление хк G <Э/?.д (—связано с нахождением оптимального вектора хк в задаче

-Ьл{-рк) = III ¡ II (У л ZA + ркх) -1л-л + ВАх < dA

ZA > О, :Г > 0.

Эта задача может оказаться неограниченной. Тогда ее двойственная является несовместной, то есть ^рк <1 dorn hа- Этот случай аналогичен недопустимости в прямой задаче (8), и тогда необходимо аппроксимировать область dorn hA с тем, чтобы дальнейшие приближения не выходили из этого множества.

Для двойственных переменных гиперплоскость, отделяющая вектор ^рк от множества dorn hа, будет задаваться опорным вектором прямых переменных хк. Этот вектор является решением задачи

mm{cAzA +ркх) Aaza + Вах < 0 0 < zA < е, 0 < х < е,

где е - вектор из единиц соответствующей размерности.

Однако, в данной ситуации возможен более простой метод построения требуемого отсечения. Неограниченность задачи (15) означает существование вектора о, что для любого допустимого вектора .г(, и скаляра в —^ +оо, х0 + ва ^ ^оо. Тогда отделение недопустимых векторов р достигается с помощью отсечения ар < 0, которое должно быть добавлено к ограничениям задачи (11) аналогично дополнительным отсечениям в прямой задаче (9).

5 Модифицированный алгоритм двойственных отсечений (МАДО)

В результате введения дополнительных отсечений допустимости МАДО для задачи (6) приобретает следующий вид:

Шаг 1 Выбрать начальные приближения х°, р°.

Шаг 2 Для начального приближения р° вычислить hß(p°). Определить начальные аппроксимации ]'ц 1 (>) = —оо, //л' (}>) = —оо. Положить к = 0, т = 0, п = 0 и ввести множества А','/' = 0, /'" = 0.

Шаг 3 ЕСЛИ текущая точка хк не оптимальна, то выполнять:

Решить задачу

ф(хк)= minJHI, (ВХ 0)

AßZß BßX ^ (Iн

x + s = хк > 0. ./• > 0

ЕСЛИ lls^ = 0, то решить задачу

Ii

„к

fB(x) = min cBzB (ВХ\

AßZß -Ь BßX fi dß h

x — X > 0. X > 0

и определить pk G dfB(xk), hB{pk), ИНАЧЕ

определить p™, /1™ = p™xk — ф(хк) и положить множество допустимости XJ = X™'1 U {х : pfx < т = т + 1

Модифицировать прямую аппроксимацию fB(x) = max{/£_1(a;), ркх - hB(pk), ркх - hB{pk)},

Решить задачу

min А/а(х) + fB(x)} (Ах)

и определить хк+1, fA(xk+1), ЕСЛИ задача

^hÄ{-pk)= min (cAzA+pkx) (АР)

Aaza + Вах < dA zA> 0, x> 0

неограничена, то определить отсечение допустимости вида х'^р < 0 и множество допустимости = Р^1 U {р : х%р < 0}, п = п + 1, ИНАЧЕ определить хк G dhA(-pk), fA(xk),

Модифицировать двойственную аппроксимацию hkA(p) = тархк - fA(xk), pxk+l - fA(xk+l)},

Решить задачу

min Ш-р) + hB(p)} (BP)

P£pä 1

u i гк

1 d

и определить pk+1, hB(pk+1), • к = к + 1, перейти к шагу 3.

ИНАЧЕ хк - решение задачи (6).

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

доказательства сходимости А ДО. Если понимать как частично определенную функцию вида

fUx) г < +оо, , е *ь

[ = +оо, х 0 Хк,

то отсечения дополнительности гарантируют и выполнение условия fß(xk) = /¿+1(a;fc) в том смысле, что если хк £ dorn fB, то хк £ Xk+i.

6 Тестовый пример

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

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

В качестве примера рассмотрена задача нахождения стационарного распределения температур в некоторой компактной двумерной области I) с известными условиями на границе Г. Функция распределения температур '/'(>. //) является решением уравнения Лапласа

д2Т д2Т

V*T=—+ —= 0, (х, y)£D (26)

с граничными условиями Дирихле на границе,

Т(х,у) = t(x,y), (х,у) £ Г. (27)

Известно, что решение Т(х,у), удовлетворяющее уравнениям (26)-(27), совпадает с функцией, минимизирующей функционал

_ 1

Х ~ 2

d

'д2т\2 (д2т\2'

дх2 ) \ ду2

da; dy, (28)

где '/'(.г. у) - функция из допустимого множества пробных функций, заданных в I). Для этой задачи пробные функции допустимы, если они непрерывны и имеют кусочно-непрерывные первые произодные. Кроме того пробные функции должны удовлетворять граничным условиям (27).

Рис. 3: Разбиение области на треугольные элементы

Классический подход применения метода конечных элементов к решению задачи (28) [3] приводит к решению линейной системы

КГ = Я, (29)

где К - матрица жесткости системы, Т - искомый узловой вектор, Я - вектор правой части системы. Детали вычисления соответствующих матриц можно найти, например, в[3].

Для достижения высокой точности решения необходимо использовать большое количество элементов, что соответственно повышает размерность матрицы жесткости системы К. Последнее время широко исследуются подходы, связанные с эффективным предобуславливанием системы линейных уравнений с учетом ее блочной структуры [2], [9], [11], [12].

В качестве тестовой решалась задача нахождения стационарного распределения температур в двумерной области, показанной на рис. 3. Использовались кусочно-линейные аппроксимации поля температуры в треугольных трехузловых лагранжевых элементах. На правой части границы данной области, в узлах 1-5, задано условие Дирихле Т0 = 50, на остальной части границы выполнено граничное условие Неймана = 0, где п - нормаль к границе. Известно, что граничные условия Неймана выполнены автоматически для функции, минимизирующей функционал (28), как естественное следствие вариационной формулировки. Точное решение задачи очевидно - температура во всех узлах области Т = Т0 = 50.

Количество узлов системы равно 62, область разбита на 100 элементов. "Портрет" матрицы жесткости системы с распределением ненулевых элементов приведен на рис. 4. Блочное представление матрицы жесткости системы очевидно. Пунктирной линией показаны границы подматриц Аа, Ва, Ва, Вв. На рис. 5 показана матрица жесткости

системы с разбиением на блоки. Типичные характеристики подзадач приведены в таблице 1.

0 10 20 30 40 50 60

Рис. 4: "Портрет" матрицы жесткости К

Блок А Блок В

Количество внутренних 20 32

переменных г

Количество связывающих 10 10

переменных х

Количество ограничений 25 37

Количество граничных 5 0

условий

Таб. 1: Характеристика подзадач тестового примера

7 Результаты вычислительного эксперимента

Модифицированный алгоритм двойственных отсечений реализован на языке С с использованием библиотеки функций линейной алгебры meschach [10]. Расчет производился на ПК типа PIII/850 под ОС Linux с использованием компилятора gcc. Текст программы (не считая meschach) составляет около 2000 операторов. Эта реализация может рассматриваться в качестве предварительной, в дальнейшем будет реализован параллелизованный вариант МАДО.

¿Л

X

0 10 19

0%-1- 0

10 10 10

09 0 |- 0

20 20

24 24 ^-%24

0 10 19 0 9

У:

10

й,

09

Ж

0 10 20 30

10

20

30

V

10 10

20 20

30 30

36 - 36 36

0 9 0 10 20 30

- 10

20

30

36

>>в

0

0

Рис. 5: "Портрет" матрицы жесткости с разбиением на блоки Аа, Ва, Вв.; Ав

Статистика решения тестового примера приведена в табл. 2, Результаты решения задачи показывают достаточно быструю сходимость по "большим" итерациям, оптимальное решение устанавливается уже после к = 5, это говорит о небольшом обмене данными между задачами (АХ) и (ВР), что способствует разработке параллелизованно-го алгоритма. На итерациях 1-5 задача ВР возвращает опорные векторы р™ отсечений допустимости, на 6-ой итерации алгоритма вектор // становится нулевым. Для решения подзадач МАДО использовался симплекс-метод с определением двойственных переменных [6]. Суммарное количество итераций симплекс-метода при решении подзадач (АХ) составило 265, в подзадачах (ВР) - 162.

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

Задача (АХ) на (к + 1)-ом шаге отличается от задачи (АХ) на к-ом шаге наличием дополнительных отсечений. Соответственно двойственная постановка задачи (АХ) на (& + 1)-ом шаге по сравнению с к-ым шагом содержит дополнительные переменные. От-

сюда следует, что оптимальное решение двойственной для (АХ) задачи на к-ом шаге является допустимым для этой же задачи на (к + 1)-ом шаге. Использование формулировки задачи (ВР) в двойственной постановке позволяет использовать ее оптимальное решение на к-ом шаге в качестве допустимого для этой же задачи на (к + 1)-ом шаге, полагая добавленные дополнительные переменные нулевыми.

Применение такой стратегии рестарта при решении подзадач позволяет уменьшить суммарное количество "малых" итераций симплекс-метода в подзадаче (АХ) до 45, в подзадаче (ВР) - до 26,

Насколько позволяет судить наш пока весьма ограниченный вычислительный опыт, алгоритм достаточно хорошо масштабируется при увеличении размеров задачи. Решались задачи и большей размерности с количеством связывающих переменных равным 50 и 102. Количество "больших" итераций сохраняется порядка числа связывающих переменных. Вместе с тем существует проблема повышения эффективности при решении подзадач, доминируют затраты на решение задач (ВХ), Возможно, что трудоемкость данной задачи может быть уменьшена при использовании проекционных алгоритмов для решения подзадач [4],[13].

к тип рк, р7 тип хк

0 ВР 0.36 0.54 0.36 0.54 0.36 -0.36 -0.54 -0.36 -0.54 -0.36 АХ 52.92 66.85 50.59 40.94 34.80 0.00 200.00 0.00 0.00 0.00

1 ВР 0.36 0.54 0.36 0.54 0.36 -0.36 -0.54 -1.00 -1.00 -1.00 АХ 78.97 58.87 45.47 37.45 32.29 178.11 82.26 0.00 0.00 0.00

2 ВР 0.36 0.54 0.36 0.54 0.36 -0.36 -0.54 -1.00 -0.54 -1.00 АХ 48.33 48.05 49.26 56.76 45.85 48.18 39.37 0.00 129.07 0.00

3 ВР 0.36 0.54 0.36 0.54 0.36 -0.37 -0.53 -0.37 -0.54 -1.00 АХ 50.32 53.42 53.50 49.99 41.13 38.16 63.65 84.07 54.92 0.00

4 ВР 0.36 0.54 0.36 0.54 0.36 -0.37 -0.52 -0.37 -0.54 -0.36 АХ 49.42 49.98 49.86 50.27 50.63 46.80 52.01 45.62 51.66 52.44

5 ВР 0.04 -0.12 -0.01 -0.08 0.17 -0.45 1.00 -0.55 1.00 -1.00 АХ 50.00 50.00 50.00 50.00 50.00 50.00 50.00 50.00 50.00 50.00

6 ВР 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 АХ 50.00 50.00 50.00 50.00 50.00 50.00 50.00 50.00 50.00 50.00

Таб. 2: Статистика решения тестового примера, к - номер итерации алгоритма двойственных отсечений, хк, рк - значения векторов связывающих и сопряженных переменных на к-ой итерации, />,"' - опорный вектор отсечения допустимости. Столбцы "тип" описывают типы решаемых подзадач

Литература

[1] Гловински Р., Лионе Ж.-Л., Тремольер Р. Численное исследование вариационных неравенств. ,\1.: Мир, 1979.

[2] Корнеев В.Г. Эффективное предобуславливание методом декомпозиции областей,-Известия вузов, Математика, 1999, No, 5, 11,

[3] Норри Д., де Фриз Ж, Введение в метод конечных элементов,-М,: Мир, 1981,

[4] Нурминский Е.А. Численные методы выпуклой оптимизации,-М,: Наука, 1991.

[5] Пшеничный Б.Н. Выпуклый анализ и экстремальные задачи,-М,: Наука, 1980.

[6] Танаев B.C. Декомпозиция и агрегирование в задачах математического программирования.-Мн.: Наука и техника, 1987.

[7] Цурков В.И. Декомпозиция в задачах большой размерности,-М,: Наука, 1981.

[8] S. Elhedhli, J.-L. Goffin, J.-P. Vial Cutting plane methods for nondifferentiable optimization. in "Encyclopedia on Optimization" C. Floudas, P. Pardalos, editors, Kluwer Academic Publishers, 1999. ftp://ftp.crt.umontreal.ca/pub/users/jlg/ElGoVi2_3.ps.gz

[9] J. Mandel Balancing domain decomposition, 1992. http://www.mgnet.org/mgnet/papers/Mandel/bdd.ps.gz

[10] Meschach library version 1.2b ftp://thrain.anu.edu.au/pub/meschach/

[11] S.V. Nepomnyaschikh, Method of splitting into subspaces for solving elliptic boundary value problems in complex-form domains. Russian ( Sov ) J. Numerical Anal, and Math. Model., 1991, Vol. 6, No. 2, 151-168.

[12] J.-P. Shao Domain decomposition algorithms, University of California, Los Angeles, PhD thesis, 1993. http://www.mgnet.org/mgnet/papers/Shao/jpshaol.ps.gz

[13] J. Gondzio, E. Sarkissian and J.-P. Vial Using an interior point method for a master problem in a decomposition approach iterative algorithms. Logilab, НЕС, Section of Management Studies, University of Geneva, Technical Report 1995.30, October 24, 1995, revised May 15, 1996. http://ecolu-info.unige.ch/ logilab/sarkissian/pointers/master.ps