О решении систем нелинейных алгебраических уравнений вариационным методом

иркутским государственный университет путей сообщения

Петрова Е. Г., Стрекаловский А. С.

УДК 519.853.4

О РЕШЕНИИ СИСТЕМ НЕЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ ВАРИАЦИОННЫМ МЕТОДОМ

1. Введение

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

Несмотря на широкий спектр методов, разработанных в этой области, проблема численного поиска решений систем нелинейных уравнений остается весьма актуальной. Так, при применении ньютоновских методов возникает трудность, заключающаяся в выборе подходящего начального приближения, обеспечивающего сходимость к решению [2, 3, 8, 11, 14]. Причем с ростом размерности системы сложность поиска стартовой точки многократно возрастает. В случае же использования вариационных методов, когда мы имеем дело с невыпуклой экстремальной задачей, для ее решения неприемлемы стандартные методы выпуклой оптимизации, имеющие свойство отыскивать в лучшем случае лишь локальные решения или стационарные точки, которые могут быть весьма далеки от глобального решения [3, 9, 12, 13].

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

В данной работе рассматривается вариационный подход к решению систем нелинейных уравнений с d.c. функциями, т.е. функциями, представимыми в виде разности двух выпуклых

функций. Класс d.c. функций ОС(Я") является достаточно широким векторным пространством, порожденным конусом выпуклых функций и замкнутым относительно часто встречающихся в оптимизации операций таких, как сумма, разность, максимум, минимум и т.д. [1, 12, 19]. О широте этого класса функций говорит то, что все дважды дифференцируемые функции, в частности, полиномы и тригонометрические функции попадают в

ОС (Я" ) . К тому же любая непрерывная функция на компакте может быть приближена с любой степенью точности функцией из ОС (в топологии равномерной сходимости) [12].

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

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

Рассматривается следующая задача:

/

где х еЯ", а / (х) - d.c. функции, т.е.

/

где gi(х),Н{(х) - выпуклые на Я" функции, I = 1,.,т. Как известно [2, 4], система уравнений может быть сведена к задаче безусловной минимизации

Ях) = ФСЛ

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

F(x) = ZU(x) |,

i=1

F (x) = ^(f, (x))2.

Нетрудно видеть, что задача (3) ввиду представления (2) оказывается в общем случае невыпуклой. Тем не менее, используя известные свойства d.c. функций [12, 19], нетрудно видеть, что, например, для функции (4) справедливо следующее d.c. представление:

где

F(x) = G(x) - H(x),

m

G(x) = 2Z max{g, (x)h, (x)}

,=1 m

H (x) = gt (x) + h (x)) -

выпуклые функции.

Итак, функция F(х), определяемая формулой (4), является d.c. функцией. Таким образом, для решения задачи

F(х) = G(x) - H(х) i min, х е Rn, (P) может быть использован аппарат d.c. минимизации [12]. Общая процедура глобального поиска минимума d.c. функции состоит из двух основных этапов: локального поиска и процедуры выхода из критической точки, основанной на условиях глобальной оптимальности (УГО) [12], с последующим повторным применением локального поиска. Так, например, для задачи (P) условия глобальной оптимальности имеют следующий вид [12]: z е Sol(P) ^ V(y,0) е Rn xR :

H(y) = 0-е, С := F(z), G( x) - <VH (y), х)>0-<VH (y), y) Vx е Rn. (8) Кроме того, УГО обладают так называемым алгоритмическим (конструктивным) свойством, позволяющим в случае их нарушения строить допустимую точку лучшую, чем исследуемая, на чем и основана процедура глобального поиска. Действительно, если при некоторых (y, 0) и х условие (8) нарушено:

G(х) <P + (VH(y), х - y),

то из выпуклости H вытекает F

или F(x) < F(z) . Так что точка x «лучше», чем z .

3. Специальный метод локального поиска

Для задачи d.c. минимизации (P) базовым элементом, «кирпичом» является, на наш взгляд, решение (линеаризованной в текущей точке xs е R") выпуклой задачи

G

x

В зависимости от выбора метода решения этой задачи («кирпича») глобальный поиск (все здание) может оказаться удачным или нет: «устойчивым» к выбору начального приближения или способным получать в лучшем случае просто допустимую точку.

Сам же локальный поиск может состоять из последовательного (как в методе простой итерации) решения задач (PLs). Зная xs, можно найти xs+1 как приближенное решение (PLs). Процесс в

этом случае сходится (по функции F = G — H) [12] к решению линеаризованной задачи

G

x

Точка x» в этом случае называется критической точкой относительно метода локального поиска.

Теорема [12]. Пусть в задаче (Р) целевая функция F ограничена снизу, а последовательность |строится по правилу

G(xs+1)-{VH(xs),xs+1)-8S < < inf {g(x) - (VH(xs), x) | x e Rn}

где Ss >0, s = 0,1,2,...,< • Тогда

s=0

i) Последовательность } удовлетворяет следующему условию:

inf G(x) - G(xs+1) + (VH(xs), xs+1 - x)}

lim

xeRn

= 0;

и) В случае сильной выпуклости функций О

4. Алгоритм глобального поиска

На основе условий глобальной оптимальности задача (Р) декомпозируется на ряд более простых задач [12]: одномерный поиск, построение аппроксимации поверхности уровня выпуклой

i=1

,=1

иркутским государственный университет путей сообщения

функции H, задающей базовую невыпуклость, проверка точек аппроксимации на пригодность к дальнейшим исследованиям и локальный поиск. В результате можно предложить следующий алгоритм глобального поиска (А1 II) (см. также [12]).

Пусть заданы начальная точка x0 , отрезок [Р_, Р+ ], где Р = inf {G(x)| x е Rn},

x

P+ = sup{G(x) | x е Rn } и числовые последовательности {тк} и {Sk}: Tk,Sk> 0, £ = 0,1,2, .. ^ ткЛо, бкЛ 0 (Аг ^ 00).

Шаг 0. Положить к := 0, xk := x0.

Шаг 1.1. Начиная с xk е Rn, посредством специального метода локального поиска построить тк -критическую точку

zk, 4 := F(zk ) <F(xk ).

Шаг 1.2. Если F(zk )< %, где % - заданная точность, то Stop, zk — приближенное решение системы уравнений.

Шаг 2. Выбрать некоторое Р е Р, Р+ ].

Шаг 3. Построить аппроксимацию

Шаг 9. Если щ (Р) < 0, то положить Р .= Р + АР е Р, Р ] и перейти на шаг 3.

Шаг 10. Если щ (Р) < 0 Vp е Р, р ] (т.е. одномерный поиск по Р завершен), то положить к .= к +1, xk+1 .= zk и вернуться на шаг 1.

Наибольшую трудность в предложенном методе вызывают шаги 3, 4 и 6, на которых возникает необходимость построить аппроксимацию поверхности уровня функции H , решить линеаризованную задачу

G

и задачу уровня (VH(v),U —v>t max, H(v) = P~Ck.

Очевидно, что в силу представления (7) линеаризованная функция оказывается негладкой, поэтому для нахождения ее минимума необходимо воспользоваться методами негладкой минимизации [5, 15].

Задача уровня (9), как будет показано далее, в некоторых случаях имеет аналитическое решение. В случае же, когда решение задачи (9) несет в себе большие вычислительные сложности (в общем случае она может быть невыпуклой), решение

д (p) = {v1 ... vNk | H(V) = p — ^k i = 1 ... N} задачи уровня и вычисление величины щ(P) за-

Шаг 4. Для каждого i = 1,..., Nk отыскать точку

—i Т\П «

u е R , являющуюся решением линеаризованной задачи

G(U) _<VH(vi), U) — Sk < inf {G(x)-<VH(V), x) | x еЯ"}.

x

Шаг 5.1. Для всех i = 1,..., N, исходя из точки U1, найти тк-критическую точку u' е Rn методом локального поиска.

Шаг 5.2. Если F(u' )< %, где % - заданная точность, то Stop, u1 — приближенное решение системы уравнений.

Шаг 6. Для каждого i = 1,..., N найти w1:

<VH (w' ), u' — w' ) + Sk> > sup{<VH(v),u' — v) | H(v) = Р — £к}.

v

Шаг 7. Положить щ(Р) := rfk(P) + Р, где

п ... .А

ГОР) := <VH(wj),u] — wJ) — G(uJ) =

А

= max{<VH (w'), u' — w') — G(u')}.

1<i< N

Шаг 8. Если щ (Р) > k, то положить

x

k+1 j

:= u и вернуться на шаг 1.

меняется сравнением значения целевой функции в полученной точке со значением .

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

Для вычисления нижней границы интервала одномерного поиска [Р_, Р+ ] решается задача выпуклой негладкой минимизации.

Поскольку функция О в данном случае не

ограничена сверху на Я", то верхняя граница Р+ вводится искусственным образом.

x

v

5. Численный эксперимент

Для тестирования описанного метода глобального поиска вначале были рассмотрены системы квадратичных уравнений следующего вида:

1г

где X е Яп, С — (и х и) -симметричные, не обязательно положительно определенные матрицы. Как известно [12], знаконеопределенную, но симметричную матрицу С можно различными способами представить в виде разности двух симметричных положительно определенных матриц С = А — В. Это влечет следующее d.c. представление квадратичных функций:

/. (х) =

< Ах, х> + (Ь, х> + (11

< В,х, х>

= Ег (х) — кг (хХ * =1,—,

откуда для целевой функции ¥(х) несложно получить разложение

¥

где G(х), Н(х) - выпуклые функции. В силу того, что для построения аппроксимации поверхности уровня функция Н должна иметь максимально простую структуру, рассмотрим следующее представление для G( х) и Н (х):

и

х) = ^ тах{< А х, х> + <Ь{, х> + di;

г=1

< Вх х> — <Ьг, х> — }, Н (х) = Ях, х),

и

где Я = ^ (А + В ) — (п х и) -симметричная по-

г=1

ложительно определенная матрица.

Отметим, что в случае системы квадратичных уравнений задача уровня (9) имеет аналитическое решение

(

и =

2

<Яиг , иг >

1

• и.

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

исходя из данных о задаче, и - число, которое находится из условия

Н(уг) = 1 <&,уг> = р — £к, г = 1,...,Кк.

В этом случаенабор точек аппроксимации зависит от текущей критической точки 2 и изменяется от итерации к итерации. К рассмотрению были предложены следующие аппроксимации:

^01 = {уг = —М,ег | г = 1,...,и},

где ег = (0,...,1,...,0) - вектора из стандартного базиса,

^ 02 = {уг = ^ 1 г = l,.•, u},

где sг - столбцы матрицы Я . Такое задание точек аппроксимации отражает структуру функции А Н (х) .

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

^03 = {уг = — ис' |г = 1,...,и},

где с, =(—1^-1,1,1,...,1).

I

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

^12 = ^01 ^ ^02 ^ ^03.

Для одномерной минимизации по / было использовано разбиение отрезка [/, ] с постоянным шагом, зависящим от размерности.

В качестве тестов были взяты системы нелинейных уравнений [6, 10, 17], решения которых заранее известны

/г =

/г = 3хг — хг—1 — 2х,+1 — кх2 + 1 х0 = хи+1 = 0 ;

3)

/х = 1 — х, /, = 10(г — 1 )(х, — х,—!)22, г = 2,..., и;

4)/ =1 — х, г = 1,3,5,.,и — 1,

/г = 10(х — х^д г = 2,4,.,и (и - чет-

ное);

5) = 3Х, ( Х,+! — 2Х, + х—0 +

(хг+1 — хг—1)2

4

х0 =0; хи+1 = 20.

иркутским государственный университет путей сообщения

Результаты численного решения этих систем, представленные в табл. 1, получены с помощью PC Intel Pentium 4 3.0 Ггц, программа написана на языке C++. Точность по значению целевого функционала бралась равной 0,001. В качестве стартовых выбиралось несколько начальных точек, в таблице представлены результаты для точек, обеспечивающих наихудшие условия для работы алгоритма. В табл. 1 и далее используются следующие обозначения: № - номер системы нелинейных уравнений; x0 - начальное приближение; n - размерность системы; F0 - начальное значение целевого функционала F(x); F - полученное значение функции (если алгоритмом получено решение с требуемой точностью, то F в таблице полагается равным нулю); St - количество итераций; PL - количество решенных линеаризованных задач; T - время работы алгоритма (минуты: секунды. доли). Для сравнения в таблице приведены результаты работы Matlab-солвера STRSCNE [18], предназначенного для решения систем нелинейных уравнений. В колонке F приведено значение целевой функции F , полученное в результате работы данного солвера, запущенного из рассматриваемой точки x0.

Отметим, что в первых двух системах со-лверу STRSCNE не удалось получить решение из предложенной стартовой точки, в то время как АГП во всех случаях получил глобальное решение.

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

v2

6) f = 2x, _ Х;_х _ xi+1 +Y (x +1, +1) ,

v = —Ц-, t = Х0 = x«+1 =0 ;

n +1

7) fi = x_1 _ 2x, + x+1 _ v  exp x,> где v= xo = xn+i =0;

n +1

8) / = x -—(У x3 + O-

Для этих систем также несложно получить d.c. разложение, учитывая, что

х 3 =1[(х 2 + (х +1)2)2 +1] -

4

1

~ [(x +1)4 + 2 x4 + (x2 +1)2],

Таблица 1

N x0 n Fo F* St PL T fstr

1 (1,-1,...,1,-1) 10 57 0 2 29 00:00.00 2,5630

20 117 0 2 46 00:05.12 2,6924

30 177 0 2 251 01:32.30 2.6962

40 237 0 2 302 03:04.30 2.6962

2 (1.....1) 10 12 0 6 72 00:00.12 0

20 21 0 10 452 00:00.64 0

30 30 0 16 532 00:31.48 7.2293

40 39 0 18 724 01:37.15 7.1053

3 (2,-2,..,2,-2) 10 7201 0 10 352 00:01.14 0

20 30400 0 18 756 00:26.18 0

30 69600 0 20 1443 07:28.42 0

40 124800 0 22 1612 10:02.34 0

4 (0,...,0) 10 5 0 38 4874 00:08.04 0

20 10 0 54 6234 00:56.46 0

30 20 0 74 10562 02:15.22 0

5 (10,..,10) 10 600 0 10 412 00:01.04 0

20 600 0 21 1205 01:58.52 0

N х0 п Ро р* РЬ Т РБТР

6 (1,-1,...,1,-1) 10 38.3 0.01229 2 167 10:00.00 0

20 78.2 0.07898 2 214 10:00.00 0

7 (1.....1) 10 38.01 0 1 9 00:01.31 0

20 78.1 0.04537 2 173 10:00.00 0

8 (2,-2,..,2,-2) 10 4.6 0 2 12 00:00.43 0

20 9 0 3 21 00:02.64 0

30 13.37 0 3 24 00:18.62 0

40 17.75 0 2 28 00:20.84 0

60 22.15 0 3 42 00:16.18 0

70 26.5 0 3 47 00:53.78 0

80 30.87 0 3 58 01:56.35 0,017

ХУ = 1(х + у)2 -х2 + у2) .

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

Уг = -И/, Н(уг) = , поэтому коэффициент и находился приближенно. В качестве точек аппроксимации выбирался набор .

В табл. 2 приведены результаты тестирования систем 6-8. Время на решение системы ограничивалось 10 минутами.

В то время как с помощью решателя STRSCNE системы 6-8 решены с необходимой степенью точности (за исключением системы размерности 80), алгоритму глобального поиска в трех случаях (см. табл. 2) не удалось найти решение с заданной точностью. Но значение целевой функции при этом достаточно уменьшилось, что позволяет использовать полученную в результате работы алгоритма точку в качестве стартовой для стандартных численных методов решения систем.

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

БИБЛИОГРАФИЯ

1. Александров А. Д. О поверхностях, предста-вимых разностью выпуклых функций // Изв. АН КазССР. Сер. Математика и механика. 1949. № 3. С. 3-20.

2. Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. Численные методы. М. : Наука, 1987.

3. Васильев Ф. П. Численные методы решения экстремальных задач. М. : Наука, 1988.

4. Вержбицкий В. М. Основы численных методов. М. : Высш. шк., 2002.

5. Демьянов В. Ф., Васильев Л. В. Недифферен-цируемая оптимизация. М. : Наука, 1981.

6. Дэннис Дж., Шнабель Р. Численные методы безусловной оптимизации и решения нелинейных уравнений. М. : Мир, 1988.

7. Ильин В. А., Поздняк Э. Г. Линейная алгебра. М. : Наука, 1978.

8. Ортега Дж., Рейнтболдт В. Итерационные методы решения нелинейных систем уравнений со многими неизвестными. М. : Мир, 1975.

9. Поляк Б. Т. Введение в оптимизацию. М. : Наука, 1983.

10. Набор тестовых систем нелинейных уравнений / А. Роозе, В. Кулла, М. Ломп, Т. Мерессоо. Таллин : Валгус, 1989.

11. Срочко В. А. Численные методы: курс лекций. Иркутск : Изд-во Иркут. ун-та, 2004.

12. Стрекаловский А. С. Элементы невыпуклой оптимизации. Новосибирск : Наука. Сиб. отд-ние, 2003.

13. Сухарев А. Г., Тимохов А. В., Федоров В. В. Курс методов оптимизации. М. : Наука, 1986.

14. Фаддеев Д. К., Фаддеева В. Н. Вычислительные методы линейной алгебры. СПб. : Лань, 2002.

иркутским государственный университет путей сообщения

15. Шор Н. З. Методы минимизации недифферен-цируемых функций и их приложения. Киев : Наукова думка, 1976.

16. Hiriart-Urruty J.-B. From Convex Optimization to Nonconvex Optimization // Nonsmooth Optimiza-tiom and Related Topics. N. Y. : Plenum Press, 1989. Pt. I : Necessary and sufficient conditions for global optimality.

17. More J. J., Garbow B. S., Hillstrom K. E. Testing unconstrained optimization software // ACM

Пономарев Д. В., Финогенко И. А.

Transactions on Mathematical Software. 1981. V. 7. P. 17-41.

18. Bellavia S., Macconi M., Morini B. STRSCNE : A Scaled Trust Region Solver for Constrained Nonlinear Equations // COAP. 2004. V. 28, N. 1. P. 31-50.

19. Tuy H. D. c. Optimization: Theory, Methods and Algorithms // Handbook of Global Optimization. Dordrecht: Kluwer. 1995.

УДК 517.9

О РАЗРЫВНЫХ СИСТЕМАХ УПРАВЛЕНИЯ С ОБОБЩЕННЫМИ ФУНКЦИЯМИ

Введение. В данной работе исследуются разрывные системы управления

Бх = / (?, х, щЦ, х),..ит (?, х)) + р(Г), (1) где / (?, х, щ,..., иот ) - непрерывная по совокупности аргументов векторная функция, скалярные кусочно-непрерывные (см. [1, с. 39]) функции щ (?, х) разрывны только на поверхностях Б = {(t,х): х) = 0}, каждая ограниченная часть которых имеет нулевую меру Лебега. Это означает, что в каждой точке (^ х) е Б существуют конечные пределы и+ = и+ (I, х) и и- = и- (^ х) функций щ х) с обеих сторон поверхности . Функция р(^) = (р1 (0,..., рп представляет собой обобщенную производную Бм>(1:) от некоторой ограниченной измеримой на каждом конечном интервале функции w(t) = п ^)).

Под решением задачи (1) понимается измеримая функция х(^), для которой существует измеримый селектор v(t) е Е (^ х^)) такой, что Бх^) = v(t) + р^), где Е (t, х) - выпуклая оболочка множества х,Ц(^х),..,ит(t,х)), и^(t,х) -отрезки с концами и+ = и+ (t, х), и- = и- (I, х) для каждых х) е Б . Если (t, х) £ Б, то Ui (t, х) состоит из одной точки щ (^ х) . Формально задачу (1) будем записывать в виде дифференциального включения

Бх е Р^, х) + р^),

понимая его решение в указанном выше смысле.

При условии p(t) = 0 получаем обыкновенное

дифференциальное включение х е Е(t, х), хорошо изученное в работах А.Ф. Филиппова [1], М.А. Айзермана, Е.С. Пятницкого [2], В.И. Уткина

[3].

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

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