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

С.К.Завриев, А.В.Федосова* МГУ им. Ломоносова, факультет ВМК

November 4, 1999

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

Рассмотрим задачу полубесконечной оптимизации

Р(У°) : найти х £ Д^,

XoPt = {х еХ° \ f(x) = min f(x)}

X° = {x£X°\g(x,y)<OVy£Y0},

X° С 7Zk, Y° С 7Zl. Функции f(x) и g(x,y) предполагаются непрерывно дифференцируемыми и выпуклыми по х на выпуклом компакте Х° для любого у £ У0, 1Z - т-мерное пространство.

В данной постановке множество У0 задает бесконечную систему ограничений, определяющих допустимое множество Х° задачи Р(У°).

Для решения задачи Р(У°) применим метод внешних аппроксимаций, который заключается в сведении исходной постановки задачи Р(У°) к решению последовательности аппроксимирующих задач P(Yn), п = 1, 2,следующего вида:

Р(У) : найти х £ Xopt[Y],

Xopt[Y] = {х£ X[Y] I f(x) = min f(x>)},

*19alina0cmc.msu.ru

X[Y] = {x £ Xй I g(x, у) < О У у £ Y}.

На основании указанного метода в [4] построен идеальный вариант алгоритма SMETH.ACTIV.sip и доказана теорема о сходимости. Рассмотрим реальный вариант этого алгоритма для конкретной выпуклой задачи полубесконечной оптимизации Р(У°).

Вспомогательная процедура SPROC.ACTIV.sip призвана формировать множество AYn = {у*} для добавления к множеству Yn:

Yn+i := Yn U /SYn,

если вспомогательная процедура не установит, что точка х уже является решением исходной задачи (команда ОСТАНОВ основного алгоритма).

Основная процедура метода SMETH.ACTIV.sip при помощи данных, полученных SPROC.ACTIV.sip, на каждой итерации решает задачу P(Yn) с конечным числом ограничений.

Для ОСТАНОВА используется параметр imax- Если imax р&з не нарушается неравенство на Контрольном шаге Процедуры SPROC-ACTIV.sip, то управление передается в основной алгоритм SMETH.ACTIV.sip на Шаг ОСТАНОВ и алгоритм заканчивает свою работу.

Метод SMETH.ACTIV.sip

Параметры. imax < ос. Шаг 0. (Начальный шаг)

Yl := 0.

Шаг 1. Находим х - решение задачи P{Yn). Шаг 2. Выполняем процедуру SPROC.ACTIV.sip с входными параметрами х и V .

Получаем AYn либо выход на ОСТАНОВ. Шаг 3. Формируем новое множество ограничений

Yn+i := Yn U /SYn.

Шаг 4.

п := п + 1. Переход к шагу 1.

ОСТАНОВ. Из процедуры SPROC.ACTIV.sip получена точка х - решение исходной задачи Р(У°).

Процедура SPROC.ACTIV.sip

Входные параметры : х £ X. Y £ Л4/(¥°). Выходные параметры : AY £ ./ИДУ )

либо ОСТАНОВ.

Параметры процедуры : 7 > 0, е > 0. Шаг 0. (Инициализация).

i := 1.

Шаг 1. Если / > imax, тогда возврат в SMETH.ACTIV.sip

на ОСТАНОВ. Шаг 2. (Пассивный случайный поиск критического ограничения).

Найти точку у-, £ У , используя равномерное вероятностное распределение на У0. Шаг 2. ACTIV. (Активный поиск критического ограничения). Применить алгоритм локального поиска для решения задачи д(х,у) -> max (метод проекции градиента),

y€Y°

стартуя из точки у; для получения точки у* £ У0, такой что :

V*i YstaAx)i 9{х, Уг) < д(х, У*).

Шаг 3. (Контрольный шаг). Если

г тах((ж, у\), ...,д(х, у*)) < 7,

то i := i + 1 и переход к Шагу 1. Шаг 4.

:= {}, где г* = arg max д(х, у*).

2. Численное исследование алгоритма для задачи упруго-пластического кручения стержня.

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

Рассмотрим тонкий цилиндрический стержень из упруго-идеал ь-нопластического материала, подчиняющегося критерию Мизеса в прямоугольной декартовой системе координат. Пусть ось жз параллельна образующей цилиндрической поверхности стержня, и в х \ х-> лежит поперечное сечение, которое задается открытой, ограниченной односвязной областью Q в R2. Пусть в исходном состоянии стержня отсутствуют напряжения. Стержень с одним закрепленным концом подвергается воздействию возрастающей крутящей пары сил.

Деформированное состояние стержня характеризуется углом закручивания на единицу длины стержня. Обозначим через и(х\. х->) функцию напряжений. С возрастанием приложенных внешних сил (и угла закручивания ) в стержне проявляются пластические свойства материала (появляются зоны пластичности). В таких зонах поведение материала перестает быть обратимым (после снятия напряжений сохраняются остаточные деформации).

Исходная краевая задача эквивалентна следующей вариационной задаче [3]: найти функцию и(х\. х-)), удовлетворяющую во всей области О неравенству Vu < 1 (для определенности ограничение сверху для Vu положено равным единице), на границе области - условию и\\- = 0. Г = 0Q. и доставляющую минимум функционалу

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

где Hq(Q) - пространство Соболева 1-го порядка с нулевым следом на границе Г [2]. Производные понимаются в обобщенном смысле.

Заметим, что функция и(х\. х->) будет удовлетворять неравенству Vu < 1 в области упругости, а равенству Vu = 1 в области пластичности.

Сведем вариационную постановку (2) к задаче полу бесконечной оптимизации. Для этого исходное бесконечномерное пространство

min J(u)

К = {и £ Щ(П) Vn < 1 п.в. в 0}

аппроксимируем по схеме Ритца, например, на базе ортогональных рядов Фурье. Соответствующее разложение выбираем с учетом граничных условий. При этом число переменных становится конечным (так как пространство из бесконечномерного стало конечномерным), а число ограничений остается бесконечным, ибо они были операторными и должны выполняться V.7- Е Q.

Таким образом, подставив соответствующее разложение в функционал J, и проинтегрировав его по сечению О, получим выпуклую (это следует из выпуклости функционала J) целевую функцию для задачи 57 Р. которая зависит только от коэффициентов ряда Фурье.

Допустимая область задачи полубесконечной оптимизации получается из подстановки разложения функции и(х\. х->) в ограничения

Vu(#i,жг)2 - 1 < 0 У(ж1,жг) £ -

Расчеты упруго-пластического кручения проведем для двух видов сечений стержня:

1) О = (0,1) х (0,1) - квадратная область;

2) О - единичный квадрат с вырезанной правой верхней четвертью.

Для первого случая используем следующее разложение:

t t

u(xi,X2) = X! Е aijsin(7rixi)sin(7TJX2)-i=ij=i

После его подстановки в функционал (1) получим целевую функцию:

о i=i j=i 7г i=1 j=1 IJ

где A - квадратная матрица основных переменных u;j. 7 = 1..... t. j = 1..... t. которые являются коэффициентами разложения функции и(х\. х->) в ряд Фурье.

Ограничения задачи выглядят следующим образом:

t t

д{А1х) = С^2 Е aijCos(7Tixi)sin(7TJX2)TTi)2+ i=ij=i

t t

(E E aijsin(7rixi)cos(7rjx2)7rj)2 - 1 < 0

i=lj=l

УЖ1 G (0,1) и x2 £ (0,1).

Таким образом, вариационная постановка задачи упруго-пластического кручения стержня с квадратным сечением (0,1) х (0,1) была сведена к выпуклой задаче полубесконечной оптимизации

Р(О) : Найти А £ X®ptJ

Xh = {A£X°\f(A)=imnf(A%

Х° = {А £ Х° д(А, х) < 0 Уж £ О}, Х° с ntxt, О С П2.

К задаче P(Q) применим алгоритм SMETH.ACTIV.sip, который был реализован для = 1.2.10. 20.

В методе проекции градиента шаг выбирался по правилу Армихо.

Для экономии времени счета и получения наиболее точных результатов, сначала задача считалась для 25 переменных с любой начальной точки с высокой заданной точностью 7 (например, 7 = 0.05). Далее решение этой задачи выбиралось в качестве начальной точки для задачи со 100 переменными (т.е. матрица А имела размерность (10x10)). Для улучшения эффективности также учитывались полученные ранее активные точки у £ ft, чтобы в задачах с требованием большей точности сократить время счета.

Задача считалась решенной, если гтах раз подряд значение

max(g(x,yl),...,g(x,y*max))

не превышало г = 1,..../ > .,

Параметр гтах положили равным 500, поэтому -f1- = 0.0001 (для варианта с 25 переменными). В критерии останова в задаче со 100 переменными 7 уменьшалось постепенно от 50 до 0.5, поэтому для полученного решения в любой точке сечения выполняется д(А, х) < 0.001.

По полученному решению были сделаны графические рисунки зон пластичности (там, где д(А,х) = 0) и упругости (где д(А,х) < 0). Было получено значение функции и(},. ) = 0.412, а тестовое и(!,. ) = 0.413 для fi = 10 [1].

Случай = 1 и /1 = 2 соответствует чисто упругому кручению и задача решается за 1 - 2 итерации. Остальные значения характеризуются наличием зон пластических деформаций. Зона упругости имеет вид креста исходящего из углов квадрата. Чем больше значение . тем тоньше становятся области упругости, и тем больше имеем соответствующих активных ограничений. На рисунке 1 приводятся

распределения зон для = 10, черная область соответствует зоне упругих деформаций.

Рисунок 1.

Для сечения с невыпуклой областью О - квадратом (0.1) х (0.1) с вырезанной правой верхней четвертью - разобьем О на три части как показано на Рисунке 2 и зададим it(x\. х->) следующим образом

и(х) =

щ, если х £ Oi, и-)- если х £ Q-2-и\ + и2- если х £ Оз, где

t t

it\ = Е Е aijsin(TTixi)sin(2TTJX2), i=ij=i

t t

v-> = E E ajisin(2Trixi)sin(TTJX2)-i=ij=i

П3 Ox

Рисунок 2.

Заданная таким образом и{х\. х->) будет непрерывной на всей области определения.

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

.... 1 ( тЧsin - 21)) sin(l(j + 2l))\

2isin(f) (sin (f(2j - I)) sin (§(2j + I)) 15 I 2j-l 2j + l

тф (sin (f(;/-2Q) sm (f(j + 2Q)-!4 j-2l j + 2l

Vising) (sin (f(2j-Q) щ (f (2j + Q), ,

ш (-Wi- + -WVi-J J +

1 / 2bm(f) /ш (f(j-2Q) ш + 2Q)>

о 2-, ajiakl -i г о7 i о7

2 i=2fc i=i fc=i 1=1 V 15 V J -21 3 + 21 j

ттк2 (sin (f (2j - I)) sin (f(2j + I))

4jlsin(f) (sin((j-2l)) sin((j + 2l))\ akl 15k { 7=21 + JT2H J +

Trjl (sin (f(2j-Q) sin (f(2j + 0m aft X ( 2i-/ + 2j + l )) +

Xf f f f шк (sin(U2i-k)) sin(U2i + k))

гф2к кф2г

Aiksinf) (sin{l(i-2k)) sm(f(i + 2fc))

тгР (sin ((2г - к)) sin ((2г + к)) :~2~ \ 2г - А; 2г + A;

2lsm(f) (sin (f (г - 2fc)) sin ((г + 2fc)) 15 i П iT2ife J J +

1 / Tif (sin (f(x-2fc)) sin ((x + 2fc))>

2 pi/±2,-rjGlk 2 V 2A- i + 2fc ,

2jsm(f) /sin (§(2г - fc)) (§(2г + °*аы 15 [ 2l=k Ш J +

тггА; /sin (f (г - 2k)) sin (f (г + 2fc))\ T ( V-2* + г 12k ) +

4гЬш(М) /am((2x-fc)) sm ((2г + к))\ \

akl 15; i 27 + Ш J J +

ЕЕ ЕЕ o,-iOwT + ацап-\2 2 ) (г2 + 12) +

z г=1 j=2lk=2il=l V 4 ZZO

9ЕЕЕЕ т + л н-\2 К2,\ {к2+32)+

1 i=2k j=l к=1 l=2j \ 4 )

1 * I I I ( sin(f)(k2 + 4/2) sm()(4fc2 +I2)

о е е е е -г--1- -;-

хЕЕ Е Е---:-+ aijalk---;-I +

i=i j=i k=2ii=2j \ J *

1 / /ш (f(2x-fc)) gin (f(2x + fc))>

2 ft p<OW 2 V 2г - + 2i + k ,

гф2к j=fi2l кф2г l=fi2j

sin ((i - 20) sin ((i + 2l))\ ik (sin ((г - 2fc))

v 7=21 1T21 J + { i=2k +

sin (f (г + 2fc))\ /sin (f {2j - I)) sin (f {2j +

i + 2& J * \ 2j -1 2j + l J +

j7 (sin (f (2г - A;)) sin (§(2г + /sm (§(j - 2/))

°*°w 2 \ 2i=k 2iTk J * i J~2l +

sin ((i + 2l))\ jl (sin (f (г - 2fc)) sin ((г + 2k))\ -J2l-J + 2 {-i=2k---iTTk- *

sin ((2j-Q) sin (f(2j + Q)\\

v 2j-* + 2i + / + ее(4(2 + 4Л + 4(2+Л)-

i=lj = l 71 IJ.

Ограничения в этом случае представляются как:

д(А,х) = (Y Ц aijCos{Kixi)sin{2K2x2)1) + =ij=i

t t 2

(Y Ц o>ijsin(7Tixi)cos(27TJX2)27TJ) - 1 < 0, если ж G Oi; =ij=i

* * 2

g(A,x) = (Y Ц ajiCos(27Tixi)sin(7TJX2)27ri) +

=ij=i

* * 2

(H H ajisin(27rixi)cos(7rjx2)7TJ) - 1 < 0, если ж G O2; =ij=i

д(А,ж) = (Y Ц ajiCos(27Tixi)sin(7TJX2)27ri-\-i=ij=i

t t 2

H H 0>ijCOs(7TiXi)sin(27TJX2)7Ti) +

=1j=l

(H H 0>jiSin(27TiXi)cOs(7TJX2)7TJ +

i=lj=l

t t 2

H H o>ijsin(7Tixi)cos(27TJX2)27TJ) - 1 < 0, если ж G O3.

=1j=l

Подобная постановка задачи полубесконечной оптимизации соответствует Р(О) с рассматриваемой областью О, и, следовательно, для ее решения применим алгоритм SMETH.ACTIV.sip.

Замечание.

Точки, лежащие на отрезках х\ = 0.5, могут относиться к области Q\ или О3, т.к. функция и(х\. Х2) является непрерывной по всей области задания О. Аналогично для точек, где х-) = 0.5.

В связи с тем, что в данном случае область О не выпуклая, метод проекции градиента реализовался каждый раз на области Qj. г = 1,2,3, куда попадала точка у;. выброшенная равномерным вероятностным распределением (Шаг 1 SPROC.ACTIV.sip). В этом случае более нужная точка у* (Шаг 1 ACTIV из SPROC.ACTIV.-sip) попадает на ту же самую которой принадлежала и у,.

Стохастический метод внешних аппроксимаций для данного примера исследовался для размерности t = 12, т.е. 144 переменных, начиная с t = 5, постепенно уменьшая требуемую точность выполнения ограничений - с 0.1 до 0.001. Данная задача считалась с теми

тах

же остальными параметрами для ц = у, 10, 20, 40.