МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ

НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ ЯДЕРНЫЙ УНИВЕРСИТЕТ «МИФИ»

Нелинейное программирование в современных задачах оптимизации

Рекомендовано УМО «Ядерные физика и технологии» в качестве учебного пособия

для студентов высших учебных заведений

Москва 2011

УДК 519.85(075) ББК 22.18я7 Н49

Нелинейное программирование в современных задачах оптимиза-

ции: Учебное пособие / Ю.В. Бородакий, А.М. Загребаев, Н.А. Крицына, Ю.П. Кулябичев, Ю.Ю. Шумилов. – М.: НИЯУ МИФИ, 2011. – 244 с.

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

Пособие предназначено для студентов и практикантов НИЯУ МИФИ, обучающихся по специальности «Прикладная математика и информатика». Оно также будет полезно инженерам и аспирантам, работающим в области оптимизации параметров технических систем различного назначения.

Пособие подготовлено в рамках Программы создания и развития НИЯУ МИФИ.

Рецензент д-р техн. наук, проф. А.Д. Модяев

ISBN 978-5-7262-1451-1	Национальный исследовательский
	ядерный университет «МИФИ», 2011

О Г Л А В Л Е Н И Е
ПРЕДИСЛОВИЕ...........................................................................................................	6
ОСНОВНЫЕ ОБОЗНАЧЕНИЯ И ПОНЯТИЯ........................................................	7
1. МЕТОДЫ НЕЛИНЕЙНОГО ПРОГРАММИРОВАНИЯ............................	13
1.1. Постановка задачи........................................................................................	13
1.1.1. Минимизация функции одной переменной...................................	14
1.1.2. Поиск отрезка, содержащего точку минимума .............................	17
1.2. Методы одномерной минимизации............................................................	18
1.2.1. Методы нахождения глобального минимума
унимодальных функций..................................................................	18
1.2.1.1. Прямые методы минимизации............................................	18

1.2.1.2Методы минимизации, основанные на использовании

		производных функций.........................................................	20
1.2.2. Методы поиска глобального минимума
	многоэкстремальных функций .......................................................		22
1.2.3. Методы минимизации унимодальных функций............................			22
1.2.3.1.		Метод золотого сечения......................................................	22
1.2.3.2.		Метод дихотомии................................................................	28
1.2.3.3.		Метод парабол
		(метод полиномиальной аппроксимации) .........................	32
Контрольные задачи.....................................................................................			34
1.3. Минимизация функций без ограничений
(безусловная минимизация) ........................................................................			34
1.3.1.	Методы нулевого порядка...............................................................		36
1.3.1.1.		Метод покоординатного спуска.........................................	36
1.3.1.2.		Метод ортонормальных направлений
		(метод Розенброка) ..............................................................	43
1.3.1.3. Метод сопряженных направлений (метод Пауэлла).........			47
1.3.2.	Методы первого порядка.................................................................		52
1.3.2.1.		Градиентные методы...........................................................	52
1.3.2.2.		Метод сопряженных градиентов........................................	62
1.3.3.	Методы второго порядка.................................................................		65
1.3.3.1.		Метод Ньютона (метод Ньютона-Рафсона).......................	65
1.3.3.2.		Сходимость метода Ньютона .............................................	67
1.3.3.3. Метод Ньютона с регулировкой шага................................			68
1.3.4. Метод переменной метрики (метод Девидона) .............................			70
Контрольные задачи.....................................................................................			75
1.4. Минимизация функций с ограничениями..................................................			76
1.4.1.	Метод штрафных функций .............................................................		76
1.4.2. Метод Фиакко и Мак-Кормика.......................................................			79
1.4.3.	Методы возможных направлений ..................................................		79
1.4.4.	Метод проекции градиента.............................................................		84
1.4.5.	Метод проекции градиента при линейных ограничениях............		87

3

1.4.6.	Метод условного градиента............................................................	90
1.4.7.	Метод линеаризации.......................................................................	94
1.4.8. Другие методы минимизации функции с ограничениями............		96
1.4.9. Способы определения начальной точки........................................		97
Контрольные вопросы и задачи..................................................................		98

2.НЕЛИНЕЙНОЕ ПРОГРАММИРОВАНИЕ

В ПРИКЛАДНЫХ ЗАДАЧАХ ОПТИМИЗАЦИИ........................................			99
2.1. Применение нелинейного программирования в теоретико-игровых
методах исследования сложных систем.....................................................			99
2.1.1. Теоретические предпосылки решения матричных игр.................			99
2.1.2. Основы метода фон Неймана........................................................			101
2.1.3.	Алгоритм фон Неймана.................................................................		103
2.1.4.	Математическое программирования
	в теории биматричных игр............................................................		105
2.1.4.1.		Биматричные игры
		Основные теоретические сведения.................................	105
2.1.4.2.		Нахождение ситуации равновесия
		в биматричных играх.........................................................	110
2.1.4.3. Метод Лемке – Хоусона. Теоретические основы............			118
2.1.4.4.		Алгоритм Лемке – Хоусона решения
		биматричных игр...............................................................	125
2.2. Оптимальное распределение нагрузки в системе ядерных реакторов...			131
2.2.1 Физическая постановка задачи.......................................................			131
2.2.2. Математическая постановка и решение
	оптимизационной задачи...............................................................		133
2.2.3. Максимально возможный эффект оптимизации.........................			143
2.3. Формирование банковского портфеля максимальной доходности........			145
2.3.1.	Основные характеристики ценных бумаг....................................		145

2.3.2.Постановка задачи формирования портфеля максимальной доходности при фиксированной величине

	риска................................................................................................	150
2.3.3. Решение задачи формирования оптимального портфеля
	с использованием множителей Лагранжа....................................	152
2.3.4. Пример задачи формирование оптимального портфеля.............		157
2.4. Определение начальных условий движения космического аппарата....		160
2.4.1. Математическая модель движения космического аппарата.......		160
2.4.2.	Выбор алгоритма решения............................................................	165
2.5. Использование метода линеаризации при решении задач перехвата
средств воздушного нападения.................................................................		172
2.5.1. Особенности полета истребителя. Профиль полета ...................		173
2.5.2.	Модель движения истребителя.....................................................	175
2.5.3. Задача определения минимального времени перехвата.............		180
2.5.4.	Выбор начального приближения..................................................	182
2.5.5. Тестовые примеры задач перехвата
	по минимальному углублению цели............................................	184
2.6. Использование методов нелинейного программирования при оценке
параметров формирующего фильтра........................................................		187
2.6.1.	Основные свойства формирующих фильтров.............................	187

4

2.6.2.	Корреляционная функция формирующего фильтра
	второго порядка .............................................................................	191
2.6.3.	Задача идентификации коэффициентов формирующего
	фильтра как задача нелинейного программирования.................	194

2.6.3.1.Алгоритм решения задачи методом Ньютона – Гаусса..195

2.6.3.2.Алгоритм решения задачи

		методом покоординатного спуска....................................	199
	2.6.4.	Пример расчета коэффициентов формирующего фильтра ........	201
2.7. Оптимальный перехват цели с одним и двумя разворотами..................			203
	2.7.1. Характеристики вертикального профиля истребителя
		перехватчика ..................................................................................	203
	2.7.2.	Характеристики горизонтального профиля
		полета истребителя........................................................................	203
	2.7.3.	Решение задач перехвата с одним разворотом
		методом штрафных функций........................................................	207
	2.7.4.	Перехват цели с двумя разворотами ............................................	209
	2.7.5.	Решение задачи перехвата с двумя разворотами
		методом штрафных функций........................................................	210
2.8.	Оптимальное размещение формуляров объектов
	на электронной карте.................................................................................		217
2.9	Оптимизация режима работы ядерного реактора в
	переменном суточном графике нагрузки с
	учетом возможности утилизации энергии. ..............................................		223
	2.9.1	Постановка задачи. ........................................................................	224
	2.9.2	Анализ оптимального режима......................................................	227
2.10	Оптимизационные задачи при наличии негерметичных
	тепловыделяющих сборок в РБМК..........................................................		229
	2.10.1	Задача выбора оптимальной очередности
		извлечения негерметичных ТВС при ограничении
		на предельно-допустимый уровень выброса активности...........	230
	2.10.2	Задача выбора ТВС для выгрузки
		по негерметичности с учетом штрафа..........................................	236
	2.10.3	Задача о выборе оптимального времени выгрузки
		негерметичной ТВС с учетом штрафа .........................................	238
СПИСОК ЛИТЕРАТУРЫ..........................................................................................			241

5

ПРЕДИСЛОВИЕ

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

Данная книга является логическим продолжением предыдущего учебного пособия «Линейное программирование в современных задачах оптимизации» и представляет собой результат многолетнего преподавания курса «Математическое программирование» на кафедре «Математическое обеспечение систем» в НИЯУ МИФИ.

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

Авторы выражают искреннюю благодарность дипломникам и аспирантам кафедры «Математическое обеспечение систем»: Бушину С.В., Зубкову Д.А., Лихачевой Е.В., Вахромееву П.В, Болицевич Ю.А., осуществивших программную реализацию прикладных задач, приведенных в главе 2.

6

ОСНОВНЫЕ ОБОЗНАЧЕНИЯ И ПОНЯТИЯ

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

1. Пусть x и y – n -мерные векторы

x			y
	1			1		,
x =			, y =
xn			yn

x т

= (x , …, x

n

)

– транспонированный вектор,

1

n

тогда < x, y >= x т y = ∑xi yi – скалярное произведение векторов x

и y .

i=1

a

… a

11

1n

2.

A =

… … …

– матрица размерности m ×n .

am1

… amn

3.

f (x) = f (x1 , …, x1n )

– функция n переменных.

Функция f (x)

имеет 1-ю и 2-ю непрерывные производные на

множестве D в n-мерном пространстве Rn .

Пусть некоторая точка

x 0 D . Тогда

f (x) может быть пред-

ставлена в виде ряда Тейлора относительно этой точки:

n

∂f

(xi − xi0 ) +

f (x1 , …, xn ) = f (x10 , …, xn 0 ) + ∑

∂x

i=1

x0

1

n

n

∂

2

f

(xi − xi0 )(x j − x0j ) + 0(

x − x 0

2 ) .

+

∑ ∑

2 i=1 j=1 ∂xi ∂x j

x

0

7

0

∂f

∂f

∂f

т

4. f (x

,

, …,

– градиент функции в точке

∂x

∂x

∂x

) =

2

1

n

x

0

x 0 D .

Градиент – это вектор, который направлен в сторону наискорейшего роста функции и ортогонален в точке х0 к поверхности

(или линии) уровня	f (x) = const (рис. O.1).
x2	линии
	уровня

max

f(x)

x1

Рис. O.1. Линии уровня f (x) = const для функции двух переменных

		∂2	f		…	∂2	f
		∂x ∂x				∂x ∂x
	1		1			1	n
5. H (x 0 ) =	…				…	…			– вещественная симметричная
		∂2	f		…	∂2	f
	∂x ∂x					∂x ∂x
		n	1			n	n

матрица (матрица Гессе).

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

Можно записать

	f (x) = f (x 0 )+ < f (x 0 ), (x − x 0 ) > +	1	< (x − x 0 ),
		2

H (x0 )(x − x0 ) > +0( x − x0 2 ) .

8

6. Пусть

	a		x		+ a			x	2
	11 1					12
a21x1 + a22 x2
	a		x		+ a		i2	x	2
		i1 1
			x	+ a				x
a	m1					m2			2
			1

+…+ a1 j x j +…+ a1n xn = b1;

+…+ a2 j x j +…+ a2n xn = b2 ;

…

+…+ aij x j +…+ ain xn = bi ;

…

+…+ anj x j +…+ amn xn = bm .

– система линейных уравнений в скалярном виде, тогда Ax = b – система уравнений в матричной форме, где A – матрица коэффи-

циентов aij (i =1, …, m; j =1, …, n), b – m -мерный вектор. Система линейных уравнений называется совместной, если она

имеет хотя бы одно решение.

Пусть r – ранг матрицы и представляет собой наибольший порядок отличного от нуля определителя матрицы A (число линейно независимых столбцов матрицы A ).

r ( A) = r ( Ab ) – необходимое и достаточное условие совместности системы линейных уравнений, т.е. ранг матрицы A должен быть равен рангу распределенной матрицы системы.

Ранг не превосходит числа неизвестных n , r ≤ n . Тогда:

если r = n , то решение системы единственно;

если r < n , то система имеет бесчисленное множество решений.

7. a0 + a1 x1 +…+ an xn ≥ 0 – линейное неравенство. Совокупность точек пространства, координаты которых удовле-

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

8. Общая задача математического программирования формулируется так:

найти максимум (минимум) функции f (x) при ограничениях

gi (x1 , …, x j , …, xn ) ≤ bi , i =			1, m	;	(O.1)
x j ≥ 0 , j =		.			(O.2)
	1, n
9

x′j = xj − xj _ min

При этом:

функция f (x) называется целевой функцией;

система неравенств (O.1) и условия неотрицательности переменных (O.2) называются системой ограничений задачи;

всякое решение задачи с учетом системы ограничений, т.е. совокупность значений переменных x1 , …, xn , удовлетворяющих ус-

ловиям (O.1) и (О.2), называется допустимым решением; совокупность точек n -мерного пространства, удовлетворяющих

системе ограничений, образует так называемую допустимую область;

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

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

Замечания. Постановка задачи математического программирования является достаточно общей. Если какая-то практическая задача, на первый взгляд, кажется иной, то, как правило, с помощью простых математических преобразований её можно свести к задаче рассматриваемого вида.

1. Если требуется найти минимум функции f (x) , то это эквивалентно

поиску максимума функции − f (x ) , т.е. min f (x) = −max [− f (x)] .
x	x

2. Если заданны неравенства вида

g(x) ≥ b ,

то простой переменой знака можно прийти к виду (О.1), а именно

−g(x) ≤ −b .

3.Если на переменные x j не наложено условие неотрицательности, а заданно ограничение

xj ≥ xj _ min ,

то, вводя замену переменных

,

получаем для новой переменной условие (О.1), т.е. x′j ≥ 0 .

10

4. Функции f (x) может иметь несколько экстремумов, а именно:

-локальный и глобальный экстремумы (пусть это будет максимум);

-функция f (x) , определенная на области D, достигает на ней гло-

бального максимума в точке x 0 D , если неравенство f (x) ≤ f (x 0 )

справедливо для любой точки x D ;
- функция f (x) достигает локального максимума в точке	x 0 D , ес-
ли неравенство f (x) ≤ f (x 0 ) справедливо для любой точки x	из некото-

рой окрестности точки x0 .

5. В математическом анализе для нахождения экстремума функций используются производные. Эти методы называются классическими методами оптимизации.

Различают задачи безусловного и условного экстремумов.

Задача на безусловный экстремум – задача максимизации функции f (x) при отсутствии ограничений. Необходимые условия экстремума в

этом случае записываются в виде системы уравнений:

∂f (x1 , …, xn ) = 0 , j =1, n .

∂xj

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

Задача на условный экстремум предусматривает ограничения в виде равенств: т.е. минимизировать f (x) при ограничениях g(x) = b.

В векторной форме g и b – m-мерные векторы.

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

m

L(x, λ) = f (x) + ∑λi [bi − gi (x)],

i=1

где λi – множители Лагранжа.

Функция L(x, λ) зависит от т+п переменных, на которые наложены

ограничения.

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

∂L

∂f

m

∂g

= 0, i =1, n ,

−∑λi

i

= 0 ;

∂x

∂x

∂x

j

j

=

j

i 1

∂L

= 0, j =1, m ,

bi − gi (x) = 0 .

∂λ

i

11

6. Классические методы оптимизации применяют лишь для решения сравнительно простых задач математического программирования, так как они имеют ряд недостатков.

• Для использования методов надо, чтобы функции f (x) и g(x)

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

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

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

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

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

12

Г л а в а 1

МЕТОДЫ НЕЛИНЕЙНОГО ПРОГРАММИРОВАНИЯ

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

Общая задача нелинейного программирования формулируется так. Найти min f (x1 , …, xn ) при ограничениях

gi (x1, …, xn ) ≤ bi , i =1, m .

Пример. Задача о размещении. Эта задача напоминает

транспортную задачу.
Пусть имеются			n пунктов потребления некоторой продукции,
причем b j ( j =		)	– объем потребления в j-м пункте; а также m
	1, n

пунктов производства. Будем считать, что для каждого i-го пункта производства известна зависимость стоимости производства fi от объема производства xi , т.е. функции fi (xi ) ,

i =1, m – это, как правило, нелинейные функции (рис. 1.1).

fi (xi )

xi

Рис. 1.1. Типичный вид зависимости стоимости производства от объема производства

Из рисунка следует, что чем больше объем, тем меньше себестоимость единицы продукции.

13

Наконец, задана матрица транспортных расходов cij ,

элементами которой являются стоимости перевозок единицы продукции из i-го пункта производства в j-й пункт потребления.

Требуется найти такие объемы перевозок xij из i-го в j-й пункт и

n

такие объемы производства xi = ∑xi j , которые обеспечивают

j=1

потребности по всем продуктам в j -м пункте назначения

m

( b j = ∑xi j ) и минимизируют суммарные расходы.

i=1

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

найти

			m n			m
min f (xi j ) = ∑∑ci j				xi j + ∑ fi (xi )
			i=1 j=1			i=1
при ограничениях
m
				j =1, m;
∑xi j = b j ,
i=1
			n
				i =1, n;
xi		= ∑xi j ,
			j=1
x	i	j	≥ 0.

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

1.1.1. Минимизация функции одной переменной

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

14

На первый взгляд кажется, что задача минимизации функции одной переменной является, довольно, элементарной. В самом деле, если функция f(x), которую нужно минимизировать на отрезке [a, b], дифференцируема, то достаточно найти нули производной, присоединить к ним концы отрезка, выделить из этих точек локальные минимумы и, наконец, среди последних найти ту точку, в которой достигается абсолютный (глобальный) минимум. Этот метод является классическим методом. Он основан на дифференциальном исчислении и довольно подробно описан в литературе.

Однако для широкого класса функций эта задача не так уж проста, и классический метод имеет весьма ограниченное применение, поскольку задача решения уравнения f ′(x) = 0 может

оказаться весьма сложной. К тому же в практических задачах часто неизвестно, является ли f (x) дифференцируемой функцией. Кроме

того, во многих практических задачах часто невозможно найти явную зависимость f (x) . Поэтому существенное значение

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

Существование локальных минимумов функции f (x) , почти

всегда затрудняет поиск глобального минимума. Поэтому многие приближенные методы минимизации применимы только тогда, когда любой локальный минимум функции f (x) является

одновременно и глобальным. Это дает гарантию сходимости методов. Если же такие сведения о функции не известны, то методы применять можно, но без гарантии сходимости.

Одним из классов функций, удовлетворяющих указанному условию, является класс унимодальных функций.

Определение. Функция f (x) называется унимодальной на

отрезке [a, b], если она непрерывна на [a, b] и существуют такие числа α и β (a ≤ α ≤ β ≤ b) , что:

1)на отрезке [a, α] при a < α, функция монотонно убывает;

2)на отрезке [β, b] при β < b, функция монотонно возрастает;

3) f ( x) = f * = min f ( x) при x [α, β], т.е. данная функция

[ a ,b ]

имеет минимум.

15

На рис. 1.1.1 приведены некоторые варианты унимодальных функций.

f (x)

f (x)

a	α	β	b	x	a	α = β	b x
f (x)					f (x)

a

α = β = b

x

a = α

β

b x

Рис. 1.1.1. Варианты унимодальных функций

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

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

16

1.1.2. Поиск отрезка, содержащего точку минимума

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

Рис. 1.1.2. Поиск минимума функции f(x)

Алгоритм.

1.Положить k =1.

2.Выбрать точку x0 и определить направление убывания функции f(x).

Для этого положить шаг h > 0 и вычислить значение функции

f(x0 + h).

Если f(x0 Если f(x0 Если f(x0 Если f(x0

+h) < f(x0), то перейти к п. 3, положив x1 = x0 + h.

+h) ≥ f(x0), то положить h = –h и вычислить f(x0 + h).

+h) < f(x0), то перейти к п. 3, положив x1 = x0 + h.

+h) ≥ f(x0), то положить h = h/2 и повторить п. 2.

3.Удвоить шаг, т.е. положить h = 2h и вычислить xk+1 = xk + h.

4.Вычислить f(xk+1).

Если f(xk+1) < f(xk), то положить k = k + 1 и перейти к п. 3.

Если f(xk+1) ≥ f(xk), то поиск прекратить и в качестве отрезка, содержащего точку минимума, выбрать отрезок [xk-1, xk+1].

17

1.2. Методы одномерной минимизации

Рассмотрим с общих позиций ряд методов, позволяющих найти минимум функции f (x) при ограничениях x [a,b] .

1.2.1. Методы нахождения глобального минимума унимодальных функций

1.2.1.1. Прямые методы минимизации

Данные методы основаны на вычислении значений функции f(x) в некоторых точках; они не используют значений производных оптимизируемой функции.

Метод перебора – простейший, но редко используемый в серьезных задачах.

Согласно этому методу отрезок [a,b] делится на п равных частей

		b − a
	точками				, i =1, n . Вычисляются значения функции
		xi = a + i	n

в этих точках, и путем сравнения определяется точка минимума хт:

	f (x m ) = min f (xi ).
	0≤i≤n
В качестве точки экстремума полагается:		x* ≈ x ,		f * ≈ f (x ) .
			m	m
При этом очевидно,	что погрешность ε в	определении точки
*	b − a
минимума x равна отрезку деления, т.е. ε =		n	.

Метод перебора, предполагающий предварительный выбор точек xi, называется также пассивной стратегией поиска точки минимума x*. На практике точки xi выбираются заранее, когда удобно провести (n +1) независимый эксперимент по измерению

значений функции f(x), а последовательное измерение этих значений трудоемко или невозможно по каким-либо причинам.

Однако использование информации о функции f(x) для выбора очередной точки xi измерения (вычисления) функции f(x), уже

полученной в предыдущих экспериментах, приводит к более эффективному поиску точки x*.

18

Методы минимизации, в которых точки xi определяются в

процессе поиска с помощью найденных ранее значений f(x), называются последовательными методами. Практически все рассматриваемые ниже методы относятся к этому классу.

Метод золотого сечения. Метод состоит в том, что исходный отрезок [a, b] уменьшается по определенному закону, постепенно

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

Рис. 1.2. Схема разбиения отрезка для поиска экстремума

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

Метод Фибоначчи. Этот метод почти полностью совпадает с методом золотого сечения, но есть два отличия:

1)отрезок делится с помощью чисел Фибоначчи –

γ0 = γ1; γn = γn−1 + γn−2 , n ≥ 2,

врезультате получаем следующую последовательность чисел:

1, 1, 2, 3, 5, 8, 13, 21,…;

2) требуется до начала работы метода задать число шагов n (так

как на первой итерации отрезок делится пропорционально γn−2 и

γn

γn−1 , а величина п изменяется в обратную сторону к нулю).

γn

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

19

Следует отметить, что как в методе «золотого сечения», так и в методе с использованием чисел Фибоначчи легко аналитически рассчитать количество вычислений функции на отрезке [a, b] , если

желаемая точность вычисления x* – ε, и, наоборот, какая будет точность при п вычислениях функции.

Метод золотого сечения имеет несколько меньшую скорость сходимости, чем метод Фибоначчи.

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

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

Идея метода такова: если на отрезке [a, b] с внутренней точкой минимума есть основание полагать, что функция f (x) достаточно

хорошо аппроксимируется многочленом (2-й, 3-й степени), то за приближенное значение x* целесообразно взять точку минимума этого многочлена.

1.2.1.2. Методы минимизации, основанные на использовании производных функции

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

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

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

«золотого сечения» и с использованием чисел Фибоначчи.

Метод касательных. Данный метод используется только для выпуклых функций и имеет простой геометрический смысл: находят абсциссу c точки пересечения касательных к графику функции f (x) , проведенных в граничных точках отрезка (рис. 1.3),

для этого нужны производные в этих точках.

20

Рис. 1.3. Определение точки с

Далее анализируется знак производной в этой точке с. При этом возможны следующие ситуации (рис. 1.4).

Рис. 1.4. Возможные наклоны касательной в точке с

Если f ′(c) > 0 , то в качестве нового отрезка берется отрезок [c, b] , если f ′(c) < 0 , то в качестве нового отрезка берется отрезок

[a, c] .

Таким образом, отрезок уменьшается от итерации к итерации. Метод Ньютона. Данный метод использует не только первую, но и вторую производные функции. При определенных условиях он обеспечивает значительно более высокую скорость сходимости к точке минимума, чем рассмотренные выше методы минимизации. Для реализации этого метода функция f(x) должна быть

выпуклой, дважды дифференцируемой функцией.

Прежде всего, выбирается начальное приближение х0 и строится последовательность:

		f	′		)
xn = xn−1	−		(x			, n =1, 2,…
				n−1
			′′	(xn−1 )
		f
				21

Вычисления заканчиваются, например, если		′	(x)		≤ ε . При
		f
неудачном выборе х0 метод может расходиться.
1.2.2. Методы поиска глобального минимума
многоэкстремальных функций

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

Метод ломаных. Этот метод может быть использован для поиска глобального минимума функции, удовлетворяющей условию Липшица на отрезке [a, b] . Функция f(x) удовлетворяет

условию Липшица на отрезке [a, b] , если существует число L > 0

(константа Липшица), такое, что

′

′′

≤ L

′

− x

′′

′

′′

[a,b] .

f (x ) − f (x )

x

для всех x , x

Геометрический смысл метода ломаных состоит в построении последовательности ломаных, приближающихся к графику функции f(x) снизу и имеющих угловые коэффициенты всех звеньев ± L .

Метод ломаных невозможно реализовать без знания константы Липшица L. Ее оценку получить можно, но иногда это представляет значительные трудности.

1.2.3. Методы минимизации унимодальных функций

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

1.2.3.1. Метод золотого сечения

Как известно [2; 14], «золотым сечением» называется деление отрезка на две неравные части, такие, что выполняется соотношение:

22

весь отрезок = большая часть . большая часть меньшая часть

Нетрудно проверить, что золотое сечение отрезка [a, b] производится двумя точками:

x1 = a + r1 (b − a) ,

x2 = a + r2 (b − a) ,

где

r1 = 3 −2 5 = 0,381966...,

r =	5 −1	= 0,618033...
2	2

На рис. 1.5 изображен пример деления отрезка [a, b] в пропорции «золотого сечения».

Рис. 1.5. Пример деления отрезка [a, b] в пропорциях «золотого сечения»

Точки х1 и х2 расположены симметрично относительно середины отрезка, при этом выполняются соотношения:

r + r =1,		r	= r2 .
1	2	1	2

Замечательно здесь то, что точка х1, в свою очередь, производит «золотое сечение» отрезка [a, x2 ] :

x2	− a	=	x1 − a		.
x − a
			x	− x
1			2	1

Аналогично, точка x2 производит золотое сечение отрезка

[x1 , b] .

Опираясь на это свойство «золотого сечения», можно предложить следующий метод минимизации унимодальной функции f(x) на отрезке [a, b] .

23

Идея метода такова. На отрезке [a, b] возьмем точки x1 и x2 ,

производящие золотое сечение отрезка, и вычислим значения функций f(x1) и f(x2) в этих точках (рис. 1.6).

Рис. 1.6. Выбор отрезка, содержащего точку минимума

Если f (x1 ) ≤ f (x2 ) , то при x > x2 функция не может убывать, т.е. точка минимума лежит на отрезке [a, x2 ] . Если же f (x1 ) > f (x2 ) , то, наоборот, минимум функции лежит на отрезке

[x1 , b] .

Итак, отрезок, на котором ищется минимум, уменьшился. Процесс повторяется сначала. При этом весьма важно то, что внутри отрезка [a, x2 ] (это справедливо для рассматриваемого

примера) содержится точка х1 с вычисленным значением f(x1), которая производит «золотое сечение» отрезка [a, x2 ].

Следовательно, на данном шаге достаточно вычислить лишь одно значение f(x) во второй, новой точке, производящей золотое сечение.

Аналогично в случае, если выбран отрезок [x1 , b] .

Алгоритм метода золотого сечения.

1. Вычислить длину отрезка [a, b] − = b − a . Взять внутри отрезка две точки:

24

xл = a + r1 , xп = a + (1 − r1 ) .

2.Вычислить f (xл ) и f (xп ) .

3.Если f (xл ) ≤ f (xп ) перейти к шагу 4, если f (xл ) > f (xп ) – к шагу 5.

4. Положить a = a ,		b = xп . Вычислить	= b − a . Если	≤ ε , то
процесс прекратить. В противном случае положить:
xп = xл ,	f (xп ) – известно;
xл = a + r1 ,	f (xл ) – вычислить.
Перейти к шагу 3.		a = xл . Вычислить
5. Положить b = b ,			= b − a . Если	≤ ε , то
процесс прекратить. В противном случае положить:
xл = xп ,	f (xл ) – известно;
xп = a +(1−r1 )	,	f (xп ) – вычислить.
Перейти к шагу 3.
Блок-схема алгоритма изображена на рис. 1.7.
Замечания. 1. После каждого шага длина отрезка уменьшается в 1/ r2
раз. После п итераций длина отрезка составляет			= r2n (b −a) .

2. Численная реализация метода обладает существенным недостатком. Так, приближенные значения чисел r1 и r2 порождают накопление ошибок, которое довольно быстро может привести к нарушению пропорций отрезков в смысле отношения «золотого сечения» и, как следствие, к нарушению условий вложенности отрезков и, тем самым, к расходимости процесса. Таким образом, метод не является устойчивым по отношению к ошибкам в определении параметров r1 и r2.

Рекомендации:

с максимально возможной для используемой ЭВМ точностью определять эти параметры и точки деления отрезка;

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

В качестве примеров использования метода «золотого сечения» рассмотрим следующие задачи.

Задача 1. Сколько шагов n метода «золотого сечения» обеспечивают заданную точность ε .

25

Вход

Ввод исходных данных: a , b , ε ,

= b − a, xл = a + r1 ,

xп = a + (1 − r1)

Вычисление fп = f (xп),

fл = f (xл)

Нет

Да

fл ≤ fп

a = xл ,

b = xп ,

= b − a,

= b − a ,

xл = xп ,

xп = x л ,

xп = a + (1 − r1 )

x л = a + r1

fл = fп,

fп = fл,

fп = f (xп)

fл = f (xл)

Нет

≤ ε

Да

x = xл, fл ≤ fп;

xп, fл > fп.

Выход

Рис. 1.7. Блок-схема алгоритма метода «золотого сечения»

26

Решение. Как отмечалось выше, условием достижения точности является выполнение неравенства

n ≤ ε .

Поскольку

n = r2n (b − a) ,

то отсюда

r n (b − a) ≤ ε,

2

и, как следствие:

ε

r n

≤

,

2

b − a

n

ε

ln(r2

)≤ ln

,

b − a

n ln r2

ε

≤ ln

,

b − a

так как r2 <1, то ln r2 < 0 . Следовательно,

ε

ε

n ≥ ln

ln r2

≈ −2,1 ln

.

− a

b

b

− a

Так как

ε

< 0 ,

то n ≈ 2 .

ln

b − a

Задача 2. Найти методом золотого сечения точку минимума функции f (x) = x2 − x на отрезке [0, 2] с точностью ε = 0,1 .

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

2x* −1 = 0 = −0, 25 .

Разрешив последнее уравнение относительно x* , получим: x = 0,5, f = −0,25.

Решим задачу с помощью метода «золотого сечения». В соответствии с рассматриваемым методом r1 = 0,382 , r2 = 0,618 .

27

Учитывая это, результаты пошаговых расчетов можно свести в табл. 1.1.

										Таблица 1.1
n		n	bn		an		xл	xп	f (xл)	f (xп)
1		2	0		2		0,764	1,236	–0,180	0,26
2		1,236	0		1,236		0,472	0,764	–0,241	–0,18
3		0,764	0		0,764		0,292	0,472	–0,207	–0,24
4		0,472	0,292		0,764		0,472	0,584	–0,249	–0,2
5		0,292	0,292		0,584		0,404	0,472	–0,241	–0,2
6		0,18	0,404		0,584		0,472	0,515	–0,249216	–0,2
7		0,112	0,472		0,584		0,515	0,541	–0,249775	–0,2
Как		следует из		таблицы,			x ≈ 0,515,		f ≈ −0,249775 , что

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

Рис. 1.8. Графическое изображение функции (задача 2)

1.2.3.2.Метод дихотомии

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

28

В основу метода положен тот очевидный факт, что производная унимодальной функции меняет знак на отрезке поиска только один раз – в точке минимума (рис. 1.9).

Рис. 1.9. Характер производной унимодальной функции

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

содержит точку экстремума.
Это – общая	идея	метода,	которая фактически сводит
рассматриваемый	метод	к поиску	′
			корня уравнения f (x)= 0 на

отрезке [a, b] известным методом дихотомии. Эта идея работает, когда минимум функции f(x) является внутренней точкой отрезка [a, b]. Если x лежит на границе, то ситуация несколько иная,

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

Метод дихотомии дает более быстрое уменьшение отрезка, чем метод «золотого сечения». После n итераций длина отрезка

составляет 0,5n (b − a) , тогда как при использовании метода

29

«золотого сечения» длина отрезка будет – r n (b − a),	r ≈ 0,61 .
2	2

Приведем таблицу коэффициентов уменьшения исходного отрезка
[a, b] после n итераций для методов дихотомии						и «золотого
сечения» (табл. 1.2).
						Таблица 1.2
n	1	2	3			10
r n	0,61	0,38	0,23	...		0,0080
2
0,5n	0,5	0,25	0,125			0,001

Основной недостаток метода заключается в вычислительных трудностях, связанных с определением производных функции f (x) .

Алгоритм метода дихотомии (поиск минимума).
1.	Задан отрезок [a, b]. Вычислить длину отрезка	= b − a .
2.	Вычислить среднюю точку отрезка xc = a + 0,5	и значение

производной

dfdx (xc ).

Если производная равна нулю (с точностью ε1 ), то вычисления прекращаются – x ≈ xc .

Если производная меньше нуля, перейти к шагу 3. Если производная больше нуля – к шагу 4.

3. Положить a = xc . Вычислить = b − a .

Если ≤ ε2 , то закончить вычисления, положить x ≈ xc . В противном случае – перейти к шагу 2.

4. Положить b = xc . Вычислить = b − a .

Если ≤ ε2 , то закончить вычисления, положить x ≈ xc .

В противном случае – перейти к шагу 2. Блок-схема алгоритма приведена на рис. 1.10.

30

Вход

Исходные данные a, b, ε1, ε2, f (x)

= b −a

xc = a + 0,5

Да

Вычисление

f ′(xc )

f ′(xc )

≤ ε1

Нет

Да

b = xc

f ′ < 0

a = xc

= b −a

Да

Нет

≤ε2

a = xc

Рис. 1.10. Блок-схема алгоритма метода дихотомии

Выход

31

1.2.3.3. Метод парабол (метод полиномиальной аппроксимации)

Идея метода весьма проста: если на отрезке [a, b] , внутри

которого лежит точка минимума x , функция f (x) достаточно хорошо аппроксимируется многочленом второй степени, то за

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

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

достаточно		малой	длины,	содержащий x (например, после
нескольких шагов по методу золотого сечения).
Пусть	известны		три	значения x1 < x2 < x3 ,		таких, что
f (x )> f (x		)< f (x	). В	этом случае x [x , x	3	]. Построим
1	2	3		1

многочлен второго порядка, который проходит через три данные точки (рис. 1.11).

Рис. 1.11. Графическое представление функции f (x) = α0 + α1x + α2x2

32

Коэффициенты α0 , α1 , α2 определяются из системы уравнений:

f (x )= α + α x + α x2 ,

1 0 1 1 2 1

f (x2 )= α0 + α1 x2 + α2 x22 ,f (x3 )= α0 + α1 x3 + α2 x32 .

Минимум полинома определяется из условия

df (x)		= α + 2α	2	xˆ* = 0 .
dx		1
		x=xˆ*

Откуда получаем

xˆ = − α1 . 2α2

Найдя α1 и α2 из линейной системы уравнений, получим

		1			1	f (x1 )
	xˆ* = −				1	f (x2 )
		2
					1	f (x3 )

x12				1	x1	f (x1 )
x22				1	x2	f (x2 )		=
x32				1	x3	f (x3 )			(*)

=1 (x22 − x32 ) f (x1 )+ (x32 − x12 ) f (x2 )+(x12 − x22 ) f (x3 ).

2 (x2 − x3 ) f (x1 )+ (x3 − x1 ) f (x2 )+(x1 − x2 ) f (x3 )

Алгоритм метода парабол.								x1 < x2 < x3 ,
1.	Найти		тройку				чисел		таких,	что
f (x1 )> f (x2 )< f (x3 )						(это можно		сделать, например, используя
алгоритм поиска отрезка, содержащего точку минимума).
2.	Вычислить оценку xˆ						по формуле (*).
3.	Если		xˆ* − x		≤ ε , закончить			процесс, положив x* = xˆ* . В
			2
противном случае – вычислить f (xˆ* ).
4.	Из чисел x ,			x ,		x , xˆ*	выбрать необходимую тройку чисел и
		1		2		3

перейти к шагу 2.

33

Контрольные задачи

Задача 1. Методом удвоения шага найти отрезок, содержащий точку минимума функции f (x) = (x −5)2 . В качестве начальной точки взять: 1) x0 = 3 ; 2) x0 = 7 . Составить блок-схему алгоритма.

Задача 2. Методом золотого сечения найти точку минимума функции

f (x) = (x −5)2	с погрешностью ε = 0,1.
Задача 3.	Методом дихотомии найти точку минимума функции
f (x) = (x −5)2	с погрешностью ε = 0,1. Сравнить скорость сходимости
метода дихотомии и «метода золотого» сечения для данного примера.
Задача 4.	Найти точку минимума функции	f (x) = (x −5)2 методом
парабол с погрешностью менее 10 %.
Задача 5.	Найти точку минимума функции	f (x) = (x −5)2 методами:

Фибоначчи, Ньютона и методом касательных с погрешностью менее 5 %. Сравнить полученные результаты.

1.3. Минимизация функций без ограничений (безусловная минимизация)

Перейдем к изучению методов нахождения минимума функции многих переменных f (x1 , …, xn ) по всему пространству R n , т.е. при условии, что при поиске экстремума каждая из переменных xi может принимать значения от –∞ до +∞. Допустим, что f (x)

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

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

Все методы безусловной минимизации сводятся к нахождению последовательности точек x0 , x1 , …, x n , значения функции в которых убывают:

34

		f (x0 ) > f (x1 ) >	> f (xn ) .		(1.1)
Эти	методы	называются	методами	спуска	(или

релаксационными методами).

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

В задаче без ограничений можно в качестве x0 взять любую точку. При этом, конечно, нужно использовать любую имеющуюся информацию о поведении функции f (x) с тем, чтобы выбрать x0

не слишком далеко от точки минимума.

После того, как какая-то точка выбрана, необходимо:

1)выбрать направление, вдоль которого предполагается расположить следующую точку;

2)решить, какой величины шаг должен быть сделан в выбранном направлении.

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

xk +1 = xk + λk Sk ,

(1.2)

где единичный вектор Sk задает направление, а λk – величину шага. Изменяя процедуру выбора Sk и λk, можно получать различные методы спуска.

Для задач без ограничений любое направление является возможным, однако не все направления приемлемы. Приемлемым называется такое направление Sk, для которого выполняется условие

f (xk ), Sk < 0 ,

(1.3)

Действительно, для достаточно малых λk можно написать разложение функции в ряд Тейлора

f (xk + λSk ) ≈ f (xk ) + λ f (xk ), Sk .

(1.4)

если этот член < 0, то f ( xk +λSk )< f ( xk )

Все методы спуска можно разбить на две группы:

1) прямые методы, или методы нулевого порядка – методы,

использующие только значения функций;

35

2)методы первого порядка – методы, использующие, кроме того, первые производные;

3)методы второго порядка – методы, использующие, кроме первых производных, еще и вторые производные.

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

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

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

1.3.1.Методы нулевого порядка

1.3.1.1.Метод покоординатного спуска

Наиболее простым способом определения направления спуска

является выбор в качестве ˆ одного из координатных векторов

Sk

e1 , e2 , …, en .

Через ei = (0, …, 0, 1, 0, …, 0) (i =1, n) обозначим единичный
i	n

вектор, у которого i-я компонента равна 1, а остальные – нулю. Иначе говоря, в методе покоординатного спуска на каждой итерации поиск точки с меньшим значением функции осуществляется изменением одной компоненты вектора x при неизменных остальных.

Существует несколько вариантов покоординатного спуска. Мы рассмотрим некоторые из них.

36

Покоординатный спуск с удвоением шага (первый вариант)

Опишем первый цикл метода, состоящий из п итераций.

Пусть		заданы точка x0 и	шаг λ0 . В точке x0	выбирают
начальное		направление S1 = e1	и величину шага λ1	способом
удвоения. Этот способ состоит в следующем:
1)	выбирают произвольное начальное значение шага λ1;
2)	если	f (x0 + λ1S1 ) < f (x0 ) ,	то полагают, λ1 = 2λ1	и процесс

удвоения шага продолжают до тех пор, пока убывание функции не прекратится;

3) если f (x0 + λ1S1 ) ≥ f (x0 ) , то выбирают λ1 = λ1 / 2 и переходят к п. 2.

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

Если в направлении e1 функция f (x) убывает, то фиксируют λ1 и переходят к следующей итерации. В противном случае выбирают направление S1 = −e1 и снова определяют величину шага λ1

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

Если в направлении − e1 функция убывает, то фиксируют

найденную величину шага λ1 и переходят к следующей итерации. На следующей итерации выбирают направление S2 = e2 и

полагают начальное значение шага λ2 = λ1, а начальную точку x1= x0 + λ1S1 (если поиск неудачный, то x1 = x0 ) и повторяют

процесс, как на первой итерации.

Цикл заканчивается при k = n , т.е. после того, как пройдет поиск по всем п направлениям S1 = ±e1, S2 = ±e2 , …, S n = ±en .

Если поиск по всем n направлениям оказался неудачным, то процесс прекращается: x* = xn .

В противном случае начинается новый цикл: Sn+1 = e1 и т.д.

37

Покоординатный спуск с удвоением шага (второй вариант)

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

Рассмотрим i-ю итерацию цикла (1 ≤ i ≤ n ), которая состоит в следующем:

1) если

	f (x + λei ) < f (x) , x – текущая точка,		(1.5)
	то полагают x = x + λei	и переходят к следующей итерации;
	2) если	f (x + λei ) ≥ f (x) ,	(1.6)
	то вычисляют
		f (x − λei ) ;	(1.7)
	3) если
		f (x − λei ) < f (x) ,	(1.8)
	то полагают x = x − λei	и переходят к следующей итерации;
	4) если	f (x − λei ) ≥ f (x) ,	(1.9)

то x не меняют (полагают x = x ) и переходят к следующей итерации.

Таким образом, осуществляется n итераций цикла.

Если неравенства (1.6) и (1.9) имеют место для всех 1 ≤ i ≤ n , то уменьшают величину λ, полагая, как правило, λ = λ / 2 , и переходят к следующему циклу, то есть повторяют все процедуры предыдущего цикла, но уже с шагом, уменьшенным в два раза.

Если неравенства (1.6) или (1.9) имеют место не для всех i, то шаг λ не дробится.

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

38

Покоординатный спуск с одномерной минимизацией

Цикл данного метода опять содержит n итераций. Каждая итерация состоит в поиске минимума функции вдоль одной из координат при неизменных остальных. Поиск минимума осуществляется одним из методов одномерной минимизации. Изобразим алгоритм в виде блок-схемы, представленной на рис. 1.12.

Вход

Ввод исходных данных

εi , ε2 , n, x0 , f (x0 ), k = 0

i=1

			Поиск с точностьюε1 компоненты xik :
k = k +1								i = i +1
			min{f (x1k +1 ,..., xik−+11 , xik +1 , xik++11 , ..., xnk +1 )}.
			xi

да i ≤ n

нет

нет

f (xk +1 ) − f (xk ) ≤ ε

да

xk = xk +1

Выход

Рис. 1.12. Блок-схема алгоритма метода покоординатного спуска

39

Недостаток метода покоординатного спуска (всех его вариантов) заключается в том, что он может «застревать», т.е. он может остановиться вдали от точки минимума и не обеспечить дальнейшего улучшения. Такая ситуация может возникнуть, если поверхности уровня целевой функции обладают острыми углами или очень изогнуты. На рис. 1.13, приведенном ниже и отображающем расположение линий уровня f (x) =C , представлен

пример такой ситуации.

Рис. 1.13. Линии уровня в методе покоординатного спуска

Если процесс решения привел в точку x0 , то каким бы малым ни брать шаг в направлении x1 или x2 , нельзя получить

уменьшение значения функции. Метод «застревает» в точке x0 . Рассмотрим несколько примеров применения метода

покоординатного спуска.

Задача 1. Построить блок-схему метода покоординатного спуска с удвоением шага (второй вариант).

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

k – номер цикла;

i – номер (внутри цикла) итерации движения по i-му направлению ( i =1, n );

j – счетчик неудачных шагов, когда x не меняется; x – текущее значение аргумента (без индекса).

40

Рис. 1.14. Блок-схема метода покоординатного спуска с удвоением шага

Задача 2. Найти методом покоординатного спуска с одномерной минимизацией по каждой координате минимум функции

f (x) = x12 + x22 − 4x1 +5

с точностью ε = 0,01.

41

Решение. Схему решения задачи можно представить в виде следующей последовательности действий.

Выберем начальную точку с координатами: x10 =1, x20 =1.

Значение целевой функции, соответствующее этой точке f 0 = 4 .

Цикл k = 0 . Итерация i =1:

				x1 − меняется;
				x2 =1,
необходимо найти min{f (x1 ) = x12 − 4x1 + 7}.
x1
Из условия экстремума:
		∂f	= 2x − 4 = 0 → x = 2 .
		∂x1		1	1
Итак:
x	= 2;			f 1 = 3,		f 1 − f 0		=1 > ε .
1
x2	=1,

Итерация i = 2 :

x1 = 2;

x2 −меняется,

необходимо найти min{f (x2 ) = x22 +1}.

x2

Из условия экстремума:

∂∂xf2 = 2x2 = 0 → x2 = 0 .

Итак:
x	= 2;	f 2	=1,		f 2 − f 1		= 2 > ε .
1
x2 = 0,
Цикл k =1 .
Соответственно:
итерация i =1:

x1 − меняется;

x2 = 0.

42

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

{ f (x1 ) = x12 − 4x1 +5}.

Минимум этой функции достигается при x1 = 2 , а само

значение критерия

f 3 = f (x

= 2, x

= 0) =1.

Итак:

1

2

x

= 2;

f 3

=1,

f

3 − f 2

= 0 < ε .

1

*

x2 = 0,

= 0;

Ответ:

x

f * =1.

1

*

= 0.

x2

Линии уровня, соответствующие рассматриваемой функции приведены на рис. 1.15.

Рис. 1.15. Линии уровня функции

f (x) = x12 + x22 −4x1 +5

1.3.1.2. Метод ортонормальных направлений (метод Розенброка)

В этом методе в каждом цикле производится поиск вдоль n взаимно ортогональных направлений. После завершения цикла с помощью процедуры Грама – Шмидта строится новая система ортогональных направлений. Преимущество этого метода состоит в следующем. Если целевая функция имеет узкий искривленный гребень (точнее, овраг), то поиск по n взаимно ортогональным направлениям эффективен тем, что результирующее направление стремится расположиться вдоль оси оврага. Это существенно ускоряет процесс по сравнению с методом покоординатного спуска.

43

Рассмотрим процедуры, выполняемые на k -м цикле метода. Пусть в начале цикла имеется точка x0k , числа λ1 , λ2 , …, λn –

ˆk

ˆk

ˆk

– единичные

векторы,

задающие

длины шагов, S1 ,

S2

,…, Sn

направления

поиска

экстремума.

(Для

первого

цикла

x01 ,

λ1 , λ2 , …, λn

выбираются произвольно, направления, как правило,

совпадают с координатными осями.)

ˆk

делается шаг λ1 λ1 .

В направлении S1

Если

k

ˆk

)

k

то

шаг считается

успешным

и

f (x0

+λ1S1

≤ f (x0 ) ,

k

k

k

ˆk

. Величина

полученная точка занимает место x0

, т.е. x0

= x0 +λ1S1

λ1 умножается на

число

α (α > 1, обычно α = 3),

т.е. теперь

λ1 = αλ1 .

ˆk

Если

k

)

k

) , то шаг считается неуспешным, точка

f (x0

+ λ1S1

> f (x0

x0k остается без изменений, а величина λ1 умножается на число β

(β < 0, обычно β = −0,5 ), т.е. λ1 = βλ1 .
Далее та	же	самая процедура повторяется для остальных
направлений	ˆk	ˆk
	S2 ,…, Sn . В результате будет получена новая точка

x0k и совокупность шагов λ1 , λ2 , …, λn тоже новая. Затем снова

ˆk	и задаем шаг либо αλ1 , либо
возвращаемся к направлению S1

βλ1 в зависимости от того, успешным или неуспешным был шаг до

этого.

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

Второй этап цикла состоит в выборе исходных данных для следующего цикла.

В качестве начальной точки следующего цикла x0k +1 берется последняя успешная точка.

В качестве шагов λi , i =1, n, – последние получившиеся значения. После этого выбираются новые направления.

44

Строится система векторов:

	ˆk	ˆk	ˆk
q1 = Λ1S1		+ Λ2 S2 +…+ Λn Sn ;
q2	=	ˆk	ˆk	;
		Λ2 S2	+…+ Λn Sn
…			ˆk
qn	=				,
			Λn Sn

здесь Λi – алгебраическая сумма длин успешных шагов по i-му

направлению.

Если изобразить условно направления на плоскости, то получим схему, показанную на рис. 1.16.

Рис. 1.16. Графическое представление метода Розенброка на одной из итераций

Как видно из рисунка, вектор q1 соединяет точку, в которой

процесс находился в начале цикла, с точкой, в которую он попал в конце цикла.

Данная система векторов q1 ,…, qn ортогонализируется с

помощью процедуры Грама – Шмидта.

В качестве начального (первого) направления принимается q1

(так как это направление располагается вдоль предполагаемого гребня функции).

Процедура ортогонализации Грама – Шмидта. Пусть имеем n

векторов q1 , …, qn . Необходимо получить n единичных

ˆk	ˆk
ортогональных векторов S1	,…, Sn .
	45

Рис. 1.17. Блок-схема метода Розенброка

46

ˆ

yi

Будем находить эти векторы по формуле:

Si =

, где yi

–

yi

вспомогательный вектор, i =

,

yi

= yi , yi

1

– норма вектора.

1, n

2

Векторы yi вычисляются по рекуррентным формулам:

y1 = q1 ,

…

i−1

q

, y

j

yi = qi − ∑

i

y j , i =1, n;

j=1

y j , y j

ˆ

= q1,

S1

…

ˆ

ˆ

i−1

qi , S j

ˆ

Si

= qi − ∑

ˆ ˆ

S j ,

i =1, n .

j =1

S j , S j

Блок-схема алгоритма может быть представлена в виде, изображенном на рис. 1.17. На рисунке использованы следующие обозначения:

k – номер цикла;

= – номер итерации при движениях по ˆi ; i 1, n S

j – счетчик неудачных шагов; x0 – это x в начале цикла;

x– текущий аргумент.

1.3.1.3.Метод сопряженных направлений

(метод Пауэлла)

Два вектора x и y называются Q-сопряженными (или

сопряженными по отношению к матрице Q),	если	x тQy = 0 , где
Q – положительно определенная матрица	( x и	y называют

взаимно ортогональными, если x т y = 0 , таким образом, понятие сопряженности является обобщением понятия ортогональности).

47

Направления S1 , S2 , …, Sn являются сопряженными относительно матрицы Q, если

SiтQS j = 0 при i ≠ j , i, j =1, n .

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

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

На чем основана идея сопряженных направлений и в чем преимущество поиска вдоль таких направлений? Дело в том, что если имеем квадратичную целевую функцию

f (x) = a + x тb + 12 x тQx ,

где x – n-мерный вектор, Q – матрица Гессе (Q – положительно определенная матрица), то минимум такой функции может быть найден ровно за n шагов от любой начальной точки, если поиск вести вдоль системы из n Q-сопряженных направлений.

На рис. 1.18 приведен пример формирования Q-сопряженных направлений для случая квадратичной целевой функции f (x) .

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

аналитически и дает результат x* = −Q−1b ?

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

испытывают на квадратичной функции.

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

48

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

Рис. 1.18. Определение сопряженных направлений (n = 2)

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

На рис. 1.19 приведена иллюстрация этого свойства для квадратичной функции f (x) .

Из этого следует, что если x1 и x2 – два минимума, располо-

женные на двух параллельных направлениях, то минимум квадратичной функции нужно искать на направлении x1 − x2 .

Причем легко убедится, что это направление является Q- сопряженным с S. Градиент квадратичной функции:

Рис. 1.19. Расположение касательных к линиям уровня в точках x 1 , x2

49

f (x) = b +Qx .

Поскольку xi – это точки минимума на направлении S, градиент функции в этих точках ортогонален к S, т.е.

S т f (x1 ) = 0 = S т (b + Qx1 ), S т f (x 2 ) = 0 = S т (b + Qx 2 ).

Вычитая второе уравнение из первого, получим

0 = S тQ(x1 − x 2 ) ,

т.е. S и (x1 − x 2 ) – Q-сопряженные направления.

Алгоритм метода Пауэлла.

1. Выбрать в качестве начальных направлений S10 , S20 ,…, Sn0 , совпадающие с координатными осями (k = 0). Выбрать начальную точку поиска x0k .

2. Определить с помощью одномерного поиска

				min f (x0k + λ1S1k ) .
				λ1
Положить	x k	= x k	+ λ* S k .		Затем последовательно осуществить
	1	0	1	1
одномерный		поиск	вдоль		направлений				S2k ,…, Snk , определив
последовательность точек
			xik = xik+1 + λ*i Sik , i =							.
								2, n
3. Вычислить направление S = xnk							− x0k		и сформировать новые
направления для следующей итерации:
	{S1k+1 , S2k+1, …, Snk+1}={S2k , S3k , …, Snk , S}.
4. Найти	λ* из условия				min f (x k		+ λS k ) (т.е. найти минимум
					λ	n	1

вдоль нового полученного направления) и взять в качестве начальной точки для следующей итерации

x0k +1 = xnk + λ* S k ,

положив k = k +1 , и перейти к п. 2.

Блок-схема алгоритма метода Пауэлла приведена на рис. 1.20.

50

		ε, n, x0
i = 1, n	S i = e i
x	= x 0 ,	i = 1
		λ *
		i
	min f (x + λi Si )
	λ
	x = x + λ*i Si
		x* ≈x
	Si = x −
		λ*i
	min f (x + λi Si )
	λ
	x = x + λ*i Si
	Sn	= S
	Si = Si+1

Рис. 1.20. Блок-схема алгоритма метода Пауэлла

51

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

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

1.3.2.Методы первого порядка

1.3.2.1. Градиентные методы

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

f (x), S < 0 .

Поэтому направление S = − f (x) является приемлемым. Более

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

Градиентным называется метод, согласно которому точка xk+1 выбирается в соответствии с соотношением

xk +1 = x k − λk f ( x k ) .

Взависимости от способа выбора шага λk можно получить

различные варианты градиентного метода.

52

1-й вариант (метод наискорейшего спуска). Согласно этому методу λk определяется путем решения одномерной задачи

минимизации

min f ( xk − λk f ( x k )) .

λk

В отдельных случаях эта задача может быть решена точно, аналитически.

Интересная особенность метода наискорейшего спуска состоит в том, что предыдущее и последующее направления поиска ортогональны.

Двумерная иллюстрация метода представлена рис. 1.21. Действительно, точка xk+1 лежит на данном направлении в

точке его касания с линией уровня f (x) = f (xk +1 ) , поскольку она найдена из условия минимума на направлении − f (xk ) . Следующим направлением поиска будет − f (xk +1 ) . Известно, что

градиент направлен по нормали к линии уровня, т.е. перпендикулярно касательной к ней.

Рис. 1.21. Метод наискорейшего спуска

Аналитически это можно показать следующим образом. Обозначим

ϕ(λ) = f (xk − λ f (xk )) .

Поскольку шаг λ выбирается из условия минимума, то

	dϕ
		= 0 .
	dλ
		λ
		k
		53

Найдем ddϕλ по правилу дифференцирования сложной функции;

т.е. исходя из зависимости ϕ[x(λ)] значение производной определяется соотношением

dϕ

=

dϕ

dx

.

dλ

Откуда

dx dλ

dϕ

=< f (xk − λk f (xk )), − f (xk ) > =

dλ

λ

k

dϕ

dx

производная

производная

dx

dλ

= − < f (xk +1 ), f (xk ) >= 0 .

Следовательно, направления ортогональны.

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

2-й вариант градиентного метода. На практике нередко довольствуются выбором шага, обеспечивающего убывание функции в направлении антиградиента ( f (xk +1 ) < f (xk )) , вместо

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

Опишем способ выбора шага λk на k-й итерации.

1.Выбирается некоторое произвольное значение λ (одно и то же на всех итерациях) и определяется точка x = xk − λ f (xk ) .

2.Вычисляется f (x) = f (xk − λ f (xk )) .

3.Производится проверка неравенства

f (xk ) − f (x) ≥ ελ(< f (xk ), f (xk ) >) ,

|| f ( x )||2

54

где 0 < ε <1 – произвольно выбранная константа (одна и та же на всех итерациях).

Возможно использование более простого неравенства:

f(x) < f (xk ) .

4.Если неравенство выполняется, то значение λ и берется в качестве искомого: λk = λ . В противном случае производится

дробление λ (путем умножения λ на произвольное число α <1 ) до тех пор, пока неравенство не окажется справедливым.

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

Направления поиска в этом варианте метода уже не будут ортогональны (рис. 1.22).

Рис. 1.22. Градиентный метод (второй вариант)

Можно предложить несколько критериев окончания итерационного процесса:

1)по разности аргументов на последовательных шагах

||x k +1−xk ||≤ ε1 ,

2)по разности значений функции

||f (x k +1) − f (xk ) ||≤ ε2 ,

3)по значению градиента

||f (xk ) ||≤ ε3 .

55

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

Например, для функции f (x1 , x2 )

∂f (xk ) = f (x1k + δ1 , x2k ) − f (x1k , x2k ) ,

∂x1 δ1

∂f (x k ) = f (x1k , x2k + δ2 ) − f (x1k , x2k ) .

∂x2 δ2

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

Теорема. Если функция f (x) ограничена снизу, ее градиент

f (x) удовлетворяет условию Липшица

||f (x) − f ( y) ||≤ R || x − y ||

при любых x, y E n , а выбор шага λk производится указанными ранее способами, то для любой начальной точки x0 градиентный метод сходится, т.е. || f (xk ) || → 0 при k → ∞ .

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

|| xk +1 − x* ||≤ q || xk − x* ||, 0 < q <1

(такая скорость сходимости называется линейной).

Величина знаменателя прогрессии q зависит от соотношения наибольшего М и наименьшего m собственных значений матрицы вторых производных f (x) . Достаточно малым знаменатель q будет

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

с высокой скоростью. Чем меньше отношение Mm , тем ближе к

единице знаменатель q и тем медленнее сходятся градиентные методы.

56

Можно дать геометрическую интерпретацию этого факта

(рис. 1.23).

а

б

Рис. 1.23. Типы линий уровня:

а – линии уровня, близкие к окружности; б – линии уровня «овражного типа»

Когда отношение Mm близко к единице линии уровня функции

f (x) = const близки к окружности. Для таких линий уровня градиентный метод сходится быстро (рис. 1.23, а).

С уменьшением отношения Mm линии уровня становятся все

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

(рис. 1.23, б).

57

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

Это происходит из-за того, что кривая поиска для таких функций обычно быстро спускается на «дно оврага», а затем начинает медленное перемещение к экстремуму, так как градиент оказывается почти перпендикулярным к оси оврага (рис. 1.24).

Рис. 1.24. Использование градиентного метода в случае функции «овражного» типа

Для ускорения сходимости при поиске минимума «овражной» функции существуют специальные методы – они называются «овражными».

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

Пример. Рассмотрим функцию f (x) = x12 + 25x22 .

Очевидно, что линии уровня данной функции вытянуты вдоль оси х1. Если ввести замену переменных: у1 = х1, у2 = 5х2, то функция

58

f ( y) = y12 + y22 будет иметь линии уровня в виде окружностей, и в

результате градиентный метод сходится за один шаг.

В реальной задаче подобрать нужные преобразования, как правило, оказывается достаточно сложно.

Блок-схема градиентного метода (2-й вариант) изображена на рис. 1.25.

Рис. 1.25. Блок-схема градиентного метода, где x0 – аргумент в начале итерации;

		– текущее значение аргумента; f0 – функция в точке x0 ; f – функция в точке
x
	;	f0′ – производная в точке		; λ0 – начальное λ для каждой итерации
x			x
		59

В данном методе окончание итерационного процесса осуществляется по условию – || f ′(xk ) ||≤ ε.

Задача. Найти методом наискорейшего спуска минимум

функции f (x) = 2x 2 + x

2

+ x x

2

+ x + x

2

с точностью ε = 0,05 (по

1

2

1

1

норме градиента). В качестве

начальной

точки

принять

точку

x 0 = (0,0) . Значения

функции

и градиента

в

этой

точке,

соответственно, будут:

f (x 0 ) = 0 ;

∂f

= 4x1 + x2

1

+1

f (x0 ) =

∂x1

=

.

∂f

= 2x2

+ x1

1

+1

∂x2

x =x0

Решение. Очевидно, что нетрудно найти точное решение данной задачи:

x * = − 1		;	x* = −	3	;	f * = −	2	≈ −0,2857 .
1	7		2	7		7

Решим задачу с использованием градиентного метода. В результате получим следующую последовательность действий.

1-я итерация:

	0	1			1		0		0		0			λ		− λ
f (x			;	x		= x		− λ f (x									,
		) =								) =			−		=
		1									0			λ		− λ

f (x0 ) ≈1,4 > ε.

Так как	x1 = −λ и x1 = −λ,		то		f (x1 (λ)) = f (λ) = 4λ2 − 2λ .
	1	2			λ ,	обеспечивающее минимум
Таким	образом,	значение
			0		1		0		1
функции по направлению f (x				)		из точки x		, будет λ =		.
					=				4
					1

Окончательно получаем для первой итерации:

60

1

−1 4

1

x

f (x

)

= −1 4 = −0,25 .

=

−

1 4

;

2-я итерация:

1

−1 4

f

(x

)

|| f (x1 ) || ≈ 0,35 ;

=

1 4

,

−1 4

−1 4

λ −1

(2)

4

x

=

1 4

=

;

−1 4 − λ

(λ +1)

−

4

2λ2

−2λ −4

′

1

1

f (λ) =

,

f (λ) =

(4λ −2) =

0 → λ =

;

16

16

2

2

−1 8

2

x

=

f (x

) = −0,281.

,

−3 8

3-я итерация:

2

1 8

2

f (x

,

|| f (x

) || ≈ 0,18 ;

) =

1 8

−

(λ +

1)

−1 8

1 8

x (3) =

8

− λ

=

3)

;

−3 8

1 8

(λ +

−

8

f (λ) = 4λ2 − 2λ −18 , 64

f ′(λ) = 641 (8λ − 2) = 0 → λ = 14 ;

3

−5 32

3

x

,

f (x

) = −0,2852 .

=

−13 32

4-я итерация:

3

−1 32

3

f (x

|| f (x

) || ≈ 0,044

< ε .

) =

1 32

,

61

Как видим, условие точности решения задачи выполнено; окончательно имеем:

	*		3		−5 32				*
x		= x					,	f (x		) = −0,2852 ,
				≈	−13 32

что, с заданной степенью точности совпадает истинным решением задачи.

1.3.2.2. Метод сопряженных градиентов

Метод сопряженных градиентов является частным случаем метода сопряженных направлений.

Вновь предположим, что f (x) – квадратичная функция вида

f (x) = a + x тb + 12 x тQx ,

где Q – положительно определенная матрица, которая в данном

случае совпадает с матрицей Гессе.

В том случае, если f (x) не является квадратичной, то Q также

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

Существуют различные варианты метода сопряженных градиентов. Одним из них является метод Флетчера – Ривса.

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

антиградиент в направление Sk , которое является	Q -сопря-
женным с ранее найденными направлениями.
Эти направления вырабатываются следующим образом:
S0 = − f (x0 ),
Sk = − f (xk ) + βk −1Sk −1 ,	(1.10)
62

λ* −

k 1

где коэффициент βk −1 выбирается так, чтобы сделать Sk						Q -сопря-
женным с Sk −1 . Иначе говоря, из условия
0 =< Sk ,QSk−1 >= − < f (xk ),QSk−1 > +βk−1 < Sk−1,QSk−1 >.
Отсюда		< f (xk ),QSk−1		>
	βk−1 =				.	(1.11)
		< Sk−1,QSk−1	>
Точка xk	находится с помощью		одномерного			поиска в
направлении Sk	из точки xk−1 , т.е.
где	xk = xk −1 + λ*k −1Sk −1 ,					(1.12)

λ*k −1	= arg(min f (xk −1 + λk −1Sk −1 )) .
	λk −1

Можно показать, что вырабатываемые с помощью данной процедуры направления S0 , S1 ,..., Sk являются Q -сопряженными.

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

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

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

f (x) = Qx + b .

Тогда, принимая во внимание итерационное соотношение (1.12), можно записать

f (xk ) − f (xk−1 ) = Q(xk − xk−1 ) = λ*k−1QSk−1 .

(1.13)

Откуда получаем

QSk −1 = f (xk ) − f (xk −1 ) .

Подставим полученное выражение в формулу (1.11) для βk −1 :

βk −1 = < f (xk ), f (xk ) − f (xk−1 ) > ,

<Sk −1 , f (xk ) − f (xk−1 ) >

63

или, раскрывая скалярные произведения и производя несложные операции, получим

1

βk−1 = < f (xk ), f (xk ) > −< f (xk ), f (xk−1 ) > .

<Sk−1, f (xk ) > − < Sk−1, f (xk−1 ) >

2

3

Рассмотрим каждый из занумерованных членов.

2-й член. Так как xk – точка минимума на направлении Sk −1 , то

градиент в этой точке ортогонален к направлению поиска (рис. 1.26). Откуда получаем:

< Sk−1, f (xk ) >= 0 .

Рис. 1.26. Ортогональность градиента f (xk ) и

направления S k −1

3-й член. Подставим вместо Sk −1 его выражение, определяемое формулой (3.10). Тогда получим

< Sk −1, f (xk−1 ) >=< − f (xk−1 ), f (xk −1 ) > +βk −2 < Sk −2 , f (xk −1 ) > .

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

< Sk−1, f (xk−1 ) >= − < f (xk−1 ), f (xk −1 ) > .

1-й член. Выразим f (xk −1 ) в соответствии с формулой (1.10), тогда получим

< f (xk ), f (xk −1 ) >=< f (xk ),−Sk−1 +βk−2 Sk−2 >=

64

= − < f (xk ), Sk−1 > +βk −2 < f (xk ), Sk−2 > .

Разрешим (1.13) относительно f (xk ) и подставим в последнее выражение. Тогда получим:

< f (x ), f (x − ) >= β − < f (x − ) + λ* − QS − , S − >=

k k 1 k 2 k 1 k 1 k 1 k 2

= β − < f (x − ), S − > +β − λ* − < QS − , S − >.

k 2 k 1 k 2 k 2 k 1 k 1 k 2

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

< f (xk ), f (xk −1 ) >= 0 .

Это означает, что векторы f (xk ), f (xk−1 ) ортогональны. Вообще можно показать, что все градиенты f (x0 ), f (x1 ), ...,f (xk ) взаимно ортогональны.

Итак, в окончательном виде имеем

βk −1

=

< f (xk ), f (xk ) >

=

f (xk )

2

.

< f (xk−1 ), f (xk−1 )

>

f (xk −1 )

2

Последнюю формулу можно применять для неквадратичных функций, так как в нее не входит матрица Q .

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

xk −1 − x* ≤ c xk − x* 2 , c > 0 .

1.3.3.Методы второго порядка

1.3.3.1.Метод Ньютона (метод Ньютона – Рафсона)

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

65

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

Метод Ньютона [2] является прямым обобщением метода отыскания корня уравнения ϕ(x) = 0 , где ϕ(x) – функция

скалярной переменной. Разложение в ряд Тейлора до первого члена позволяет переписать уравнение в следующем виде:

0 = ϕ(xk +1 ) ≈ ϕ(xk ) + ϕ′(xk )(xk +1 − xk ) .

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

xk +1 = xk − ϕ′(xk ) .

ϕ (xk )

Нас интересует n-мерная задача оптимизации, сводящаяся фактически к определению корня уравнения f (x) = 0 . Разложение в ряд Тейлора в этом случае дает:

0 = f (xk +1 ) ≈ f (xk ) + Q(xk )(xk +1 − xk ) ,

где Q(xk ) – матрица Гессе. Отсюда

xk +1 = xk − Q−1 (xk ) f (xk )

при условии, что существует обратная матрица Q−1 (xk ) . Эта

формула и определяет метод Ньютона минимизации функции. Если функция квадратичная

f (x) = a + x тb + 12 xT Qx ,

где Q – положительно определенная матрица то, исходя из произвольной начальной точки x0 с помощью полученной выше формулы можно получить следующую точку:

x1 = x0 − Q−1 (b + Qx0 ) = −Q−1b .

Эта точка является точкой минимума квадратичной функции. Таким образом, для квадратичной функции метод Ньютона оп-

ределяет точку минимума за один шаг.

66

1.3.3.2. Сходимость метода Ньютона

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

случае функция будет увеличиваться в направлении Q−1 (xk ) f (xk ) ,

а не уменьшаться. Кроме того, матрица Гессе на протяжении всего процесса должна быть невырожденной, так как необходимо существование обратной матрицы.

Если матрица удовлетворяет этим условиям и, кроме того, является ограниченной и удовлетворяет условию Липшица, то

существует некоторая окрестность точки минимума x* , такая, что для любой начальной точки x0 из этой окрестности метод Ньютона сходится с квадратичной скоростью:

x

− x*

≤ c

x − x*

2

, c ≥ 0 .

k −1

k

Таким образом, для сходимости метода начальная точка x0 должна выбираться достаточно близко к искомой точке минимума x* .

Подводя итог, можно следующим образом сформулировать достоинства и недостатки метода.

Недостатки:

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

сходимость метода зависит от выбора начальной точки x0 . В

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

Достоинства:

вблизи точки минимума метод обеспечивает более быструю сходимость, чем градиентные методы;

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

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

67

1.3.3.3.Метод Ньютона с регулировкой шага

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

xk +1 = xk − ρk Q−1 (xk ) f (xk ) ,

где 0 < ρk ≤1 (обычному методу Ньютона соответствует ρk =1 ).

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

Существуют два способа выбора параметра ρk . Первый способ заключается в следующем:

1) выбирается ρk =1 и вычисляется точка x = xk + pk , где pk = −Q−1 (xk ) f (xk ) ;

2)вычисляется f (x) = f (xk + pk ) ;

3)производится проверка неравенства

f (x) − f (xk ) ≤ (ερk < f (xk ), pk >) , 0 < ε < 12 ,

4) если это неравенство выполняется, то в качестве искомого значения ρk берем ρk =1 . В противном случае производится

дробление ρk до тех пор, пока неравенство не выполнится.

При втором способе значение ρk выбирается из условия минимума функции в направлении движения:

f (xk − ρk Q−1 (xk ) f (xk )) = min f {xk − ρQ−1 (xk ) f (xk )} .

ρ≥0

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

квадратичной, трудоемкость каждой итерации увеличивается. Однако при первом способе выбора параметра ρk трудоемкость каждой итерации несколько увеличивается.

68

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

Замечания. Как было отмечено в п. 1.3.3.2, на сходимость метода Ньютона большое влияние оказывает матрица Гессе. Если она не является положительно определенной, то метод расходится. Если она положительно определенная (т.е. выпуклая f (x) ), то метод сходится в

окрестности точки минимума.

Определить, является ли матрица положительно определенной можно, используя либо критерий Сильвестра, либо собственные числе матрицы Гессе. Матрица будет положительно определенная, если все угловые миноры матрицы положительны (критерий Сильвестра) или все собственные числа матрицы положительны.

В качестве примера рассмотрим следующие задачи. Задача 1. Рассмотрим функцию

f (x) = 2x12 + x22 + sin(x1 + x2 ) .

Проверить является ли для данной функции матрица Гессе положительно определенной?

Решение. Градиент функции имеет вид:

f (x) = 4x1 + cos(x1 + x2 ) . 2x2 + cos(x1 + x2 )

Матрица Гессе:

Q(x) =

4 − sin(x1 + x2 )

− sin(x1 + x2 )

.

− sin(x

+ x

2

)

2 − sin(x

+ x

2

)

1

1

Угловыми минорами будут:

1 = 4 −sin(x1 + x2 )

и

2 = 8 − 6sin(x1 + x2 ) > 0 .

Так как

sin(x1 + x2 )

≤1 , то

1 > 0 и 2 > 0 .

69

Следовательно, согласно критерию Сильвестра Q(x) –

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

Задача 2. Проверить положительную определенность матрицы Гессе функции:

f (x) = 4x12 − x22 − 2x1x2 + 6x1 − x2 − 2 .

Решение. Легко установить, что матрица Гессе данной функции будет:

Q =		8	− 2		.
		− 2	− 2

По критерию Сильвестра матрица Гессе не является положительно определенной, так как её главный определитель

= −16 − 4 = −20 < 0 .

Проверим полученный результат, используя собственные числа матрицы Гессе:

	8 − λ		− 2			= λ2 − 6λ − 20 = 0 ,
	− 2 − 2 − λ
λ1 =		6 +	116	> 0 , λ2 =			6 −	116	< 0 .
			2					2

Как видим, не все собственные числа положительны. Следовательно, матрица Гессе – не положительно определенная. Метод Ньютона расходится.

1.3.4. Метод переменной метрики (метод Девидона)

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

70

умножений при обращении матрицы размера

[n ×n]

приблизительно пропорционально n3, что представляет собой большую величину уже при относительно небольших значениях n.

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

Основная итерационная формула метода имеет вид

xk +1 = xk − λk ηk f (xk ),

(1.14)

где ηk – матрица, аппроксимирующая обратную матрицу Гессе, а λk – шаг, который выбирается путем минимизации функции в

направлении (− ηk f (xk )).

Кстати, формула (1.14) является общей для градиентных методов и метода Ньютона. В методе наискорейшего спуска роль ηk играет единичная матрица, а в методе Ньютона – обратная

матрица Гессе.

В рассматриваемом методе матрица ηk

вычисляется по

рекуррентной формуле:

+ Ak − Bk ,

где

ηk +1 = ηk

(1.15)

(

)

т

η

(

)

тη

т

x

x

k

g

k

g

k

A =

k

k

,

B =

k

,

k

(

)т

k

(

)тη

x

g

k

g

k

g

k

k

k

xk = xk +1 − xk ,

gk = f (xk +1) − f (xk ).

Вкачестве начального значения η0 для организации рекуррентного процесса (1.15) принимается η0 = E , где E –

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

71

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

f (x) = a + x тb + 12 x тQx

с положительно определенной матрицей Гессе Q .

Оказывается, что роль матрицы Ak состоит в том, чтобы обеспечить сходимость матрицы ηk +1 к обратной матрице Гессе Q−1 , а роль матрицы Bk в том, чтобы обеспечить положительную определенность матрицы ηk +1 и в пределе исключить влияние произвольного задания матрицы η0 .

Действительно,

η0 = E ,

η1 = E + A0 − B0 ,

η2 = η1 + A1 − B1 = E + ( A0 + A1 ) − (B0 + B1 ),

kk

ηk +1 = E + ∑ Ai − ∑ Bi .

i =0 i =0

В случае квадратичной функции, сумма матриц Ai ( i =1, k ) при

k = n −1 равна обратной матрице Q−1 , а сумма матриц Bi равна

матрице E . Это легко доказать, если в начале показать с помощью математической индукции, что вырабатываемые методом Девидона векторы направлений x0 , x1 , …, xn−1 образуют множество Q-

сопряженных векторов.

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

n−1

Покажем, что ∑ Ai = Q−1 . Из соотношений (1.15), учитывая,

i=0

что f (x) = Qx + b , получим:

72

gk = f (xk +1 ) − f (xk ) = Q(xk +1 − xk ) = Q xk .

Принимая во внимание последнее выражение, можно записать:

n−1

n−1

x

i

( x

i

) т

∑ Ai

= ∑

.

( xi )

т

i=0

i=0

Q xi

Рассмотрим выражение

n−1

n−1

xi ( xi )т Q x j

∑ Ai

Q x j

= ∑

.

(1.16)

( xi )

т

Q xi

i=0

i=0

Так как векторы xi ,

x j

(i, j =

; i ≠ j) взаимно сопряженные,

1, n

то все слагаемые в формуле (1.16), кроме j-го слагаемого, равны нулю.

Таким образом, получаем:

n−1		x j ( x j		)т Q x j
∑ Ai Q x j =					= x j .
		( x j )	т	Q x j
i=0
				n−1
Очевидно, такое возможно, только если ∑Ai = Q−1 .
				i=0

n−1

Также нетрудно доказать, что ∑Bi = E .

i=0

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

Если сравнить метод переменной метрики с методом сопряженных градиентов, то оказывается, что первый обеспечивает существенно более быструю сходимость, чем второй, но в большей мере подвержен влиянию ошибок вычислений. Поэтому часто используют метод переменной метрики со следующей модификацией: через конечное число шагов (обычно n ) осуществляют обновление матрицы ηk , т.е. полагают ηk = E , и

процесс начинается, как бы, сначала.

На рис. 1.26 приведена блок-схема метода переменной метрики с периодическим обновлением матрицы ηk .

73

вход

Исходные данные: f ( x ), f ( x), x0 , ε

k = 1

ηf0=' =E f ( x0 )

x0 − аргумент в начале итерации; x − аргумент в конце итерации;

f0′ = f ( x0 ); f ′ = f (x ) .

λ* − значение λ, доставляющее одномерный минимум;

k − номер итерации.

:

Найти λ*

λ* = arg min f (x0

− ληf0 ' )

π

x = x0 − λ*ηf 0 '

f '

x * = x

≤ ε

выход

k = n

η = E

k = k +1

x = x − x0

g = f ( x)'− f0 ' (x0 )

Вычислить новое η

k= k +1

x0 = x f0 ' = f '

Рис. 1.26. Блок-схема метода переменной метрики

74

Контрольные задачи

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

Задача 1. Найти минимум функции

f (x) = 4(x1 − 5)2 + (x2 − 6)2 .

Поиск минимума осуществить с использованием метода Розенброка. В качестве начальной точки положить: x0 = (8, 9)T , значение ε = 0.6 для

остановки алгоритма. Составить блок – схему алгоритма. Задача 2. Найти минимум функции

f (x) = 4(x1 −5)2 + (x2 − 6)2 .

Поиск минимума осуществить с использованием метода Пауэлла. В качестве начальной точки положить: x0 = (8, 9)T , значение ε = 0.1для

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

Задача 3. Найти минимум функции

f (x) = 4(x1 −5)2 + (x2 −6)2 .

Поиск минимума осуществить с использованием метода наискорейшего спуска. В качестве начальной точки положить: x0 = (8, 9)T , значение ε = 0.1для остановки алгоритма. Составить блок – схему алгоритма. Сравнить полученные результаты с результатами решения данной задачи методом Розенброка и методом Пауэлла.

Задача 4. Найти локальный минимум функции

f (x) = 2x12 + x2 x1 + x2 2 .

Поиск минимума осуществить с использованием градиентного метода

с постоянным шагом,	число итераций M =10 . В качестве начальной
точки положить: x0	= (0.5, 1)T Сравнить полученные результаты с

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

Задача 5. Найти локальный минимум функции

f (x) = 2x12 + x2 x1 + x2 2 .

Поиск минимума осуществить с использованием метода наискорейшего спуска с оптимизацией величины шага по направлению спуска, число итераций M =10 . В качестве начальной точки положить: x0 = (0.5, 1)T Сравнить полученные результаты с результатами решения данной задачи при использовании градиентного метода с постоянным шагом.

75

1.4. Минимизация функций с ограничениями

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

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

Множество X называют допустимой областью, а точку x X –

допустимым решением.			min( f (x)) при
Итак, задача состоит в том, чтобы найти
ограничениях
gi (x) ≤ 0, i =		.	(1.17)
	1, m

1.4.1. Метод штрафных функций

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

Например, если скалярная переменная x ограничена условием x ≤1 , то можно использовать следующую замену переменной:

x = sin θ или x = cos θ.

Если переменная x ограничена условием x ≥ 0 , то снять ограничения можно заменой x = y2 . При этом переменная y не ограничена.

Если a ≤ x ≤ b , то возможна замена: x = a + (b − a) sin 2 θ. Эти приемы имеют ограниченное применение, тем не менее ими пренебрегать нельзя.

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

76

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

F(x,l) = f (x) + ∑ψi (gi (x),li ) ,

i

где l

– некоторый

векторный

параметр

l = {li },

i =1, m ;

ψi (gi (x), li ) , i =

,

– функция

штрафа, которая

обладает

1, m

следующим свойством:

при

gi (x) ≤ 0,

i =1, m;

0

lim (ψi (gi (x), li )) =

в

противном

случае.

i

+ ∞

l →∞

) в области X

При таком задании штрафных функций F(x,l

близко к

f (x) , а вне области X

эта функция принимает большие

значения.

Идея метода штрафных функций состоит в том, чтобы вместо задачи (1.17) рассматривать задачу минимизации функции F(x,l ) при больших li .

В общем случае F(x, l ) строится так, чтобы она была гладкой и

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

Например, штрафная функция может быть построена как взвешенная сумма квадратов невязок

			m
F(x,l		) =	f (x) + ∑liϕ2 (gi (x)) ,
где			i=1
					__
gi (x), если gi (x) ≥ 0, i =					1, m;
ϕ(gi (x)) =				__
0, если gi (x) < 0,				i =1, m.
			77

Замена решения задачи (1.17) минимизацией функции F(x,l )

при больших l позволяет приблизиться к решению исходной задачи.

Однако здесь возникают следующие вычислительные трудности.

1. Если функции gi (x) – невыпуклые, то F(x,l ) также не будет выпуклой по x . Поэтому она может обладать локальными минимумами. Так как все изученные нами методы предназначены для нахождения локального минимума, то при плохом начальном приближении x0 будет найден локальный минимум функции

F(x, l ) , не совпадающий с минимумом исходной задачи.

Если функции gi (x) – выпуклые, то F(x,l ) также будет

выпуклой, и данная проблема устраняется.

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

3. Функция F(x, l ) в точках x , для которых gi (x) = 0 , при некоторых l не имеет вторых производных, т.е. градиент в этих

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

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

78

1.4.2. Метод Фиакко и Мак-Кормика (метод барьерных функций или метод внутренней точки)

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

Итак, надо найти min f (x) при ограничениях gi (x) ≤ 0 , где функции gi (x) – выпуклые и допустимая область не пуста.

Составим вспомогательную функцию

		m	1
	F(x, k) =	f (x) − k∑		, k > 0 .
		i=1 gi (x)

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

величина. При k → 0 минимум функции F(x, k) стремится к минимуму функции f (x) с ограничениями gi (x) ≤ 0 .

Для повышения точности важно выбирать малые значения k . Однако при малых k небольшие изменения x приводят к резким изменениям функции F(x, k) . Другими словами, изменения

градиента F(x, k) вблизи границ допустимой области становятся

более резкими. Это затрудняет поиск минимума и значительно снижает ценность метода.

1.4.3. Методы возможных направлений

Методы возможных направлений – наиболее исследованный класс методов выпуклого программирования (задачи выпуклого программирования – это такие задачи нелинейного программирования, в которых целевые функции и функции ограничений выпуклы). Подробное описание данных методов было произведено в работах Зойтендейка [2].

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

79

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

Рассмотрим сущность этих методов.

Имеем задачу нелинейного программирования: найти min f (x) при ограничениях

gi (x) < 0, i =1, m ,

где gi (x) , f (x) – непрерывно дифференцируемые функции.

Последовательность точек будем, как всегда, строить по итерационной формуле

xk+1 = xk + λk Sk .

Как известно, имея допустимую точку x , удовлетворяющую всем ограничениям, необходимо принять два решения:

1) выбрать направление S , которое должно быть возможным и приемлемым, т.е. на этом направлении должны лежать точки, принадлежащие допустимой области X (возможность), и функция f (x) должна в этом направлении убывать (приемлемость);

2) решить, какой величины шаг должен быть сделан в выбранном направлении S .

Вообще говоря, существует много возможных и приемлемых

направлений.	S , т.е. такого, в котором
При выборе «лучшего» направления
функция убывает в наибольшей	степени, минимизируют
< f (x), S >.

Функция убывает в направлении S , если это скалярное произведение меньше нуля, т.е. направление S образует острый угол с антиградиентом (рис. 1.27).

Рис. 1.27. К выбору допустимых направлений

80

Будем полагать, что в точке x имеется хотя бы одно активное ограничивающее уравнение g(x) = 0 , при этом точка x лежит на

границе допустимой области. В противном случае нет проблем с выбором направления – это антиградиент.

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

Если ограничивающее уравнение линейно, то «лучшее» направление S получается проекцией вектора – f ( x ) на

многообразие, характеризуемое линейным ограничивающим уравнением(рис. 1.28).

− f ( x )

g(x) = 0

уменьшение f (x)

x

X − допустимая область

Рис. 1.28. Проекция антиградиента на активное ограничение

Если активное ограничивающее уравнение является нелинейным, проблема выбора направления становится сложнее. Помимо условия < f (x), S >< 0 (движение в сторону уменьшения

f (x) ), должно удовлетворяться условие < g(x), S >< 0 (движение

всторону уменьшения g(x) ), где g(x) – активное

ограничивающее неравенство (рис. 1.29).

Возможными являются направления S , для которых < g(x), S >< 0 (при условии выпуклости g(x) ).

Оказывается, что при выборе возможного направления нужно принимать во внимание не только ограничения, которые в данной точке выполняются точно, но и те, которые выполняются «почти точно», так как в противном случае может произойти сколь угодно

81

сильное измельчение шага вдали от точки минимума. При этом процесс не обязательно будет приводить в точку минимума. Эти соображения включены в алгоритм выбора направления.

g(x)

x

Касательная плоскость к g(x) = 0

g(x) = 0

S S

X (g(x))

Рис. 1.29. Выбор возможных направлений при нелинейном ограничивающем уравнении

Возможное и приемлемое направление S получается путем решения следующей вспомогательной задачи условной оптимизации: найти max{ϕ(σ, S) = σ} при ограничениях:

< gi (x), S > +σ ≤ 0, i Iε ;

<f (x), S > +σ ≤ 0 .

<S, S >=1 .

В данной задаче Iε = {i : −ε ≤ gi (x) ≤ 0; 0 < ε <<1}, т.е. множество Iε включает в себя индексы «почти активных»

ограничений, значения которых находятся в ε-близости от границы

( − εi ≤ gi (x) ≤ 0 ).

Очевидно, если полученное в результате решения вспомогательной задачи максимальное значение σ > 0 , то смещение в направлении S приводит к уменьшению значения функции f (x) , что следует из второго условия в ограничениях

вспомогательной задачи, и не нарушает никаких ограничений основной задачи, что следует из первого условия в ограничениях вспомогательной задачи. Последнее ограничение – это условие нормировки вектора S . Чаще всего это ограничение заменяют

82

следующим: −1 ≤ s j ≤1, j =1, n , s j – элемент вектора направления

S. Таким образом, вспомогательная задача становиться задачей линейного программирования, для решения которой существуют конечно-шаговые методы (симплекс-методы). Неизвестными параметрами в этой задаче линейного программирования являются σ, и элементы вектора S ( s1 , s2 , …, sn ).

Рассмотрим один из алгоритмов метода возможных направлений.

В качестве начального приближения x0 может быть выбрана любая точка множества X , а ε0 выбирают из интервала (0, 1].

Пусть в результате k-й итерации вычислены xk и εk . Опишем k +1 итерацию.

1.Решить вспомогательную задачу и вычислить σk ≥ 0 и Sk (если Iε пусто, т.е. xk – внутренняя точка, то Sk совпадает с антиградиентом).

2.Если σk ≥ εk , то перейти к вычислению величины шага λk .

Это можно сделать одним из двух способов:
а) решая	задачу		одномерной				минимизации			функции
ϕ(λ) = f (xk ) + λSk , причем				λ		должна		лежать	в	интервале
0 ≤ λ ≤ λдоп . Здесь λдоп – расстояние от точки xk									до ближайшей
по направлению		Sk	точки	границы			множества		X ,	т.е. точка
y = xk + λдопSk		принадлежит		границе			множества		X	(если луч
xk + λдопSk	не пересекается с границей, а полностью принадлежит
множеству X , то λдоп = +∞ );
б) число	λk можно выбрать так, чтобы функция f (x) была в
точке xk+1 меньше, чем в точке xk .
Например, в качестве λk				можно взять наибольшее из чисел,
удовлетворяющих соотношениям:
f (xk ) − f (xk			+ λk Sk ) ≥		1	λk σk ,		0 ≤ λk ≤ λдоп .
				2			xk+1 = xk + λk Sk ,
После определения			λk вычислить							положить
εk +1 = εk и прейти к шагу 1.
				83

3.	Если	0 < σk < εk , то положить	xk+1 = xk ,	εk+1 = γεk , где γ
удовлетворяет условию 0 < γ <1, и перейти к шагу 1.
4.	Если	σk = 0 , то вычислить	σ*k , решив	вспомогательную

задачу с ε = 0 . Если σ*k = 0 , то процесс поиска точки минимума

закончен	( x* = xk ). В	противном				случае положить	xk +1 = xk ,
εk +1 = γεk	и перейти к шагу 1.
Замечание. Условие		σ*k = 0	является необходимым и достаточным
условием	оптимальности	точки			–	это доказывается с	помощью
				x
специальной теоремы [2].

1.4.4.Метод проекции градиента

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

выпукло,	для отыскания направления в точке xk напрашивается
мысль спроектировать точку yk		= xk − νk f (xk ) ( νk		– некоторое
фиксированное положительное		число)	на множество X		и	в
качестве	направления спуска	взять	Sk = pk − xk ,	где	pk	–

проекция точки yk на допустимое множество X (рис. 1.30). После этого надо осуществлять спуск вдоль полученного направления. Выпуклость множества X гарантирует, что такое направление Sk является возможным. Проекция точки y k на множество – ближайшая к y k точка этого множества.

Рис. 1.30. Определение направления в методе проекции градиента

84

Итак, метод проекции градиента состоит в вычислении

проекции pk	точки yk = xk − νk f (xk )		на множество X и в
выборе шага	λk таким образом,	чтобы в точке		xk+1 = xk − λk Sk ,
где Sk = pk − xk , выполнялось		условие	движения к минимуму
( f (xk+1 ) < f (xk ) ).
В зависимости от способа выбора шага λk				можно получить

различные варианты метода проекции градиента (можно использовать одномерную минимизацию и т.д.).

Прекращения процесса является условие pk = xk , которое является необходимым и достаточным условием того, что точка xk

– точка минимума (доказывается с помощью специальной теоремы).

Например, pk = xk , если f ( xk ) = 0 , так как тогда yk = xk , и нет перемещения из точки xk . Аналогичная ситуация возникает,

если xk лежит на границе и градиент ортогонален границе допустимой области (рис. 1.31).

− f ( x )

Рис. 1.31. Направление градиента совпадает с направлением градиента к границе допустимой области

Для того чтобы найти проекцию точки yk на множество X , необходимо решить задачу минимизации квадратичной функции yk − x 2 на множестве X (т.е. x X ) – тем самым будет найдена ближайшая к yk точка множества X . В общем случае эта задача

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

множество X представляет собой шар или параллелепипед в E n ,

85

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

Рассмотрим пример явного определения проекции точки на множестве X . В качестве области X будем рассматривать

замкнутый шар радиуса r с центром в точке O в пространстве E n . Как известно, уравнение шара имеет вид

n

∑xi2 = r 2 ;

i=1

тогда допустимая область X описывается неравенством

n

∑xi2 ≤ r 2 .

i=1

Пусть y – некоторая точка, для которой надо найти проекцию p на поверхность шара. Очевидно, что данная задача эквивалентна задаче минимизации функции

φ(x) = y − x 2; x X ,

при этом
	p = arg min φ(x) .
	x X

Данную задачу можно записать в виде: найти

n

min ∑( yi − xi )2

i=1

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

n

∑xi2 ≤ r 2 .

i=1

Решение этой вспомогательной задачи находиться в явном виде:

	если
y,
	ry
p =
	n
	∑yi2
	i=1

n

∑yi2 ≤ r2 ;

i=1

n

, если ∑yi2 > r2 .

i=1

86

x = f (xk ) . В

1.4.5. Метод проекции градиента при линейных ограничениях

Будем полагать, что ограничивающие уравнения линейны, т.е. имеют вид:

Ax −b ≤ 0 ,

где A – матрица ( m ×n ); b – m-мерный вектор.

Назовем гранью многогранника X (допустимой области) подпространство, определяемое любой совокупностью активных ограничений.

Предположим,	что при заданном			x X	грань, содержащая x ,
определяется уравнением
		gi (x) = 0, i I ,
где I ={i1 , i2 ,…, i p } , p < m .
Предположим,	что	любая	подматрица		размерности ( p ×n)
матрицы A имеет ранг		p , и пусть AI – подматрица, составленная
из p строк	матрицы		A ,	соответствующих активным

ограничениям. Тогда можно предложить следующий алгоритм поиска минимума функции f (x) при линейных ограничениях.

Алгоритм (один шаг, после того как найдена точка xk X ).

1. Вычислить проекцию градиента на грань, содержащую xk по

формуле
	т	т	−1
x = I n	− AI	(AI AI	)	AI f (xk ) ,
						(AI AIт )−1
где In – единичная	матрица		(n ×n) .		Матрица
существует, так как по условию задачи ранг матрицы						AI равен
числу строк в ней.
Если xk лежит внутри X ,		то множество I			пусто, в квадратных

скобках остается одна единичная матрица In и

этом случае данный метод оказывается обычным градиентным методом.

87

2. Если

x ≠ 0 ,

то найти

λдоп ,

такое, что

λдоп

x определяет

максимальное перемещение в направлении

− x , которое может

быть сделано,

не выходя за пределы X ;

λдоп

может быть найдено

путем решения задачи

λдоп = max{λ:

(xk − λ

xk X }.

После этого найти λ* как

решение

задачи

одномерной

оптимизации:

min f (x

− λ

x

λ* = arg

k

k

0≤λ≤λдоп

и вычислить xk +1 = xk − λ* xk .

Если

λ* < λдоп ,

то совокупность

активных

ограничивающих

уравнений остается без изменений.

Если λ* = λдоп , то граница области X будет достигнута и новое

ограничивающее

уравнение

станет

активным. Тогда, если

g j (xk+1) =0 ,

j I , то надо включить в I

это

j , определить новую

матрицу AI

и вернуться к шагу 1.

Если

x = 0 , то необходимо вычислить вектор μ (размерности

p ) исходя из соотношения

μ

= ( A A т )−1

A f (x

k

) .

I

I

I

Если

все

компоненты

μi

≤ 0 ,

i =

,

то

xk

оптимально

1, p

( xk = x* ), и процесс вычисления прекращается.

В противном случае исключить из AI

строку, соответствующую

наибольшей положительной компоненте вектора μ ,

и вернуться к

шагу 1.

Дело в том, что компоненты вектора

μ

эквивалентны

множителям Лагранжа, и точка xk

является оптимальной точкой в

том случае, когда все μi ≤ 0 .

88

Однако, если некоторые компоненты вектора μ положительны,

то функция f (x)	может уменьшаться		в том случае если
перемещение от xk	осуществляется в направлении уменьшения
gi (x) , т.е. в направлении, где		gi (x) < 0 .	При этом ограничение
gi (x) становится неактивным, i		исключается из I .
Это можно проиллюстрировать рис.1.32:

f (x )

x
− f ( x )	x

Рис. 1.32. Особенности метода проекции градиента

Алгоритм работает так, что проекции градиентов, расположенных симметрично относительно грани, одинаковы (рис. 1.33).

f ( x )

1

g ( x ) = 0

2	−	f ( x )

Рис. 1.33. Пояснения к использованию метода проекции градиента.

Следовательно, если f ( x ) ≠ 0 , то ситуация, когда x = 0 , может возникнуть в двух случаях:

89

вслучае 1 точка x является точкой оптимума, следовательно, вычисления необходимо прекратить;

вслучае 2, как видно из рис. 1.33, антиградиент направлен внутрь области X , следовательно, его проектировать на грань не нужно.

Когда согласно алгоритму из AI будет убрана строка, соответствующая ограничению gi (x) = 0 , в формуле для проекции

останется In и получится, что

x = f ( x ) , т.е. направление

поиска будет направлено внутрь области X в сторону убывания функции.

1.4.6. Метод условного градиента

Идея метода условного градиента состоит в выборе направления спуска на основе линеаризации функции f (x) относительно

текущей точки xk . В результате линеаризации получаем линейную функцию fL (x) вида:

f L (x) = f (xk ) + f (xk )(x − xk ) .

Пусть yk – точка, обеспечивающая минимум функции fL (x) на множестве X . Тогда направление Sk поиска экстремума исходной функции f (x) на множестве X можно определить следующим образом:

Sk = yk − xk ,

при этом основная итерационная формула имеет стандартный вид:

xk +1 = xk + λk Sk .

Таким образом, для определения направления Sk необходимо решить задачу минимизации линейной функции fL (x) на

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

90

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

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

Существуют различные способы выбора шага λk в методе условного градиента. Рассмотрим некоторые из них.

1. Величина шага λk может выбираться из условий: 0 ≤ λk ≤1;

ϕ(λk ) = min f (xk + λSk ).

λ

Очевидно, что такая задача может быть решена каким-либо методом одномерной минимизации.

2.Можно задавать λk =1 и проверять условие монотонного убывания функции f (x) при переходе от точки xk к точке xk+1 , а именно: f (xk+1 ) < f (xk ) . Если это условие не выполняется, то необходимо дробить шаг λk до тех пор, пока не выполнится условие монотонности.

3.Величина шага определяется из условия λk = αi0 , где i0 –

минимальный номер среди номеров i ≥ 0 , удовлетворяющих условию

f (xk ) − f (xk + αi Sk ) ≥ αi ε η(xk ) .

Здесь α, ε – параметры метода, 0 < α <1; 0 < ε <1 . Величина η(xk ) является минимумом линейной функции fL (x) , т.е.

η(xk ) = min f L (x) = min < f (xk ), (x − xk ) > .
x X	x X
	91

Необходимо отметить, что величина η(xk ) в условии выбора параметра шага берется по модулю, так как η(xk ) ≤ 0 . Это следует из того, что при x = xk вспомогательная линейная функция f L (x) x=xk = 0 , следовательно, значение η(xk ) , обеспечивающее

минимальное значение fL (x) , меньше или равно нулю.

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

Процесс вычислений заканчивается в точке xk , если η(xk ) = 0 . Это может случиться, например, если f (xk ) = 0 . Из этого условия следует, что если

min < f (x k ), x − xk >= 0 ,
x X
то для всех x X будет выполняться
< f (x k ), x − xk > ≥ 0 .	(1.18)

Последнее, видимо, является необходимым и достаточным (для выпуклой функции достаточным) условием того, чтобы точка xk являлась минимумом функции f (x) на множестве X (тоже

выпуклом).

Это условие является обобщением условия стационарностиf (x k ) = 0 для задачи минимизации функций на множествах.

Проиллюстрируем условие < f (x k ), x − xk > ≥ 0 графически

(рис. 1.34).

x

x

xk

f (xk )

Рис. 1.34. Графическая интерпретация скалярного произведения

< f (x k ), x − xk > ≥ 0

92

Оно означает, что угол между градиентом в точке xk и любым вектором x − xk для всех точек x X должен быть острым. Однако, если точка xk лежит внутри области X , то этого не может быть. Следовательно, условие (1.18) справедливо, только в случае:

f (x k ) = 0 .

Если точка xk лежит на границе области, то минимум в этой точке достигается лишь тогда, когда антиградиент перпендикулярен границе и направлен из области X (рис. 1.35).

Рис. 1.35. Определение оптимального направления в случае, когда xk лежит на границе области

В этом случае углы между градиентом и любым вектором (x − xk ) будут острыми (рис. 1.35, а).

Если антиградиент отклоняется от данного направления, то найдутся точки x , для которых эти углы будут тупыми (рис. 1.33, б).

Геометрический смысл метода условного градиента поясняет рис. 1.36.

Рассмотрим пример использования метода условного градиента, который решается в явном виде.

Пример. Найти min f (x) при ограничении

n	~	)	2	≤ r	2	.
∑(xi − xi

i=1

93

Линииуровня функции fL =< f (xk ), x − xk > ортогональны f (xk )

	X	f (x0 )
		x0
	f (x1 )	x2
	x1
	y0	y1

Линииуровня < f (x0 ), x −x0 >=const

Рис. 1.36. Геометрический смысл метода условного градиента

Как видим, допустимая область представляет собой шар радиуса

~

r в n-мерном пространстве с центром в точке x .

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

Определить точку yk , являющуюся решением задачи. Найти

	min{ f L (x) =< f (xk ), x − xk >}
	x
	при ограничении
	n	~		)	2	≤ r	2	.
	∑(xi − xi
	i=1
	Решение оказывается следующим:
	~			f (xk )
	yk = x	− r							.
			|| f (xk ) ||

1.4.7. Метод линеаризации

Пусть требуется найти min f (x) при ограничениях

gi (x) ≤ 0, i =1, m .

94

Заменим в точке xk функцию f (x) и все ограничения gi (x) на

линейные путем линеаризации:

f L (x) =< f (xk ), x − xk >;

gi (xk )+ < g(xk ), x − xk >≤ 0, i =1, m.

В результате получим задачу линейного программирования: найти

min {f L (x) =< f (xk ), x − xk >}

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

gi (xk )+ < g(xk ), x − xk >≤ 0, i =1, m.

Можно было бы решение линеаризованной задачи взять в качестве следующего приближения, как это делается в методе Ньютона для решения систем нелинейных уравнений. К сожалению, как правило, этот путь не приводит к решению исходной задачи, так как вспомогательная задача линейного программирования часто не имеет решения. Поэтому необходимо наложить некоторые ограничения на приращение вектора x в точке xk , чтобы решение линеаризованной задачи в точке xk не

уходило слишком далеко от xk , оставаясь в такой окрестности xk ,

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

Таким образом, после того как получена точка xk ,						в качестве
точки xk+1 берется решение задачи минимизации: найти
	1			2
min ϕk (x) =			x −x k		+ βk < f (xk ), x − xk >
	2

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

gi (xk )+ < g(xk ), x − xk >≤ 0, i =1, m ,

где параметр βk > 0 .

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

95

1.4.8.Другие методы минимизации функций

сограничениями

Квадратичное программирование. Задачей квадратичного программирования называется задача, в которой квадратичная функция минимизируется на многогранном множестве: найти

				1
	min f (x) = a + b тx +				x тQx	,
				2

( Q – положительно определенная матрица) при ограничениях:

Ax ≤ G .

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

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

Метод Ньютона. Этот метод можно применять для задач с ограничениями; при этом направление поиска находится не явно, а получается в результате решения вспомогательной задачи. При этом формула итерационного процесса имеет стандартный вид:

xk +1 = xk + λk Sk ,

где вектор Sk = yk − xk и определяется из решения задачи минимизации на множестве X квадратичной функции

Ψk (x) =< f (xk ), x − xk > + 12 (x − xk )т Q(xk )(x − xk ) .

В последней формуле Q(xk ) – матрица Гессе.

Шаг λk может быть выбран различными способами, в

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

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

Метод покоординатного спуска. Описанный ранее метод покоординатного спуска не трудно модифицировать применительно

96

к задаче минимизации функции	на параллелепипеде: найти
min {f (x)} при ограничениях
ai ≤ xi ≤ bi ,	i =		.
		1, n

1.4.9.Способы определения начальной точки

Врассмотренных ранее методах минимизации требовалось в качестве начальной точки x0 выбрать некоторую допустимую

точку x0 X . Для X , таких, как, например, параллелепипед, шар, гиперплоскость, указать такую точку нетрудно. Однако нередко задача определения x0 является весьма непростой.

Например, если

X = {x :	qi (x) = 0,	i =				},	(1.19)
				1, m
то для определения точки x0 X		нужно					решать систему
уравнений (в общем случае нелинейных).
Чтобы найти какую-либо точку множества
X = {x :	qi (x) ≤ 0,	i =			},		(1.20)
			1, m

придется решать систему неравенств, что представляет собой весьма серьезную задачу. Если ограничения qi (x) ≤ 0 линейны, то

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

Задачу нахождения точки x0 , принадлежащей множествам

(1.19) или (1.20), можно переформулировать в виде задачи минимизации.

В случае множества (1.19) введем функцию

m

F(x) = ∑qi2 (x), x E n ,

i=1

ав случае множества (3.19) – функцию

97

F(x) = ∑m (max{qi2 (x), 0}) p , x E n , p ≥1,

i=1

и рассмотрим задачу минимизации: найти min F(x), x E n .

Эта задача решается любым из методов безусловной минимизации.

Если множество X не пусто, то условие x X равносильно условию F(x0 ) = min F(x) = 0 = F * .

Если F * > 0 , то X – пустое множество.

Контрольные вопросы и задачи [60]

Задача 1. Найти минимум функции	f (x) = x2 − 4x		при ограничении
g1 (x) = x −1 ≤ 0 методом штрафных функций.
Задача 2. Найти минимум функции	f (x) = x2	+ x2	при ограничениях
	1	2

g1 (x) = −x1 +1 ≤ 0, g2 (x) = x1 + x2 − 2 ≤ 0 методом штрафных функций.

Задача	3.	Найти минимум	функции	f (x) = 1 (x +1)3 + x					2	при
					3	1
		g1 (x) =1 − x1 ≤ 0,	g2 (x) = −x2	≤ 0
ограничениях					методом барьерных
функций (применить обратную штрафную функцию).
Задача 4. Для задачи минимизации функции				f (x) = (x −6)2 +(x −2)2
					1			2
при ограничениях − x1 + 2x2 ≤ 4,			3x1 + 2x2 ≤12, − x1 ≤ 0,				− x2 ≤ 0
указать множество возможных направлений в точке B = (2,3) .
Задача 5	Найти минимум функции f (x) = (x				− 4)2 + (x		2	−5)2		при
				1

ограничении g1 (x) = x1 + x2 −1 = 0 методом проекции градиента.

98

Г л а в а 2

НЕЛИНЕЙНОЕ ПРОГРАММИРОВАНИЕ

ВПРИКЛАДНЫХ ЗАДАЧАХ ОПТИМИЗАЦИИ

2.1.Применение нелинейного программирования

втеоретико-игровых методах исследования сложных систем

2.1.1. Теоретические предпосылки решения матричных игр

Матричной игрой будем называть антагонистическую игру [39], в которой каждый игрок имеет конечное множество стратегий.

Антагонистической игрой J,{Si}i J ,{Hi}i J называется

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

J {1,2},

H2(s) H1(s),

s S .

Такая игра задается прямоугольной матрицей

A

aij

,

i

1,m

,

j

1,n

,

где aij – значение выигрыша игрока 1, если он выбрал свою i-ю

стратегию, а игрок 2 выбрал свою j-ю стратегию. Данная игра называется игрой размерности (m n), а матрица А – платежной

матрицей.

Если игрок 1 выбирает номер строки i , а второй – j , то в результате выбора игроками независимых стратегий игрок 2 платит игроку 1 выигрыш aij . Следовательно, игрок 1 может гарантиро-

вать себе выигрыш не менее значения: v	maxmin a	ij	– нижняя
1	i j

цена игры, а второй гарантировать себе проигрыш не более вели-

чины v2	minmaxaij – верхняя цена игры. В общем случае
	j i

v1 v v2 .

99

Если min max aij = max min aij			= a	*	* = v , то в такой игре мини-
j i	i	j	i	j
максные стратегии i* и	j*	являются оптимальными, так как реше-

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

i* и j* называются чистыми стратегиями.

Если min max aij ≥ max min aij , то в игре нет ситуации равнове-
j i	i j

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

Смешанной			стратегией игрока 1 называется вектор
X = (x , ..., x	m	)т ,	удовлетворяющий условиям:
1

m
∑xi	=1, xi ≥ 0, i =	1, n	,
i=1

где число xi – вероятность, с которой игрок 1 выбирает свою i -ю

чистую стратегию. Величина

m n
H (X ,Y ) = ∑ ∑aij xi y j = X т AY	(2.1)

i=1 j=1

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

Ситуация (X * , Y * ) называется ситуацией равновесия в сме-

шанном расширении матричной игры, если для любых X и Y выполняются неравенства

	H (X ,Y * ) ≤ H (X * ,Y * ) ≤ H (X * ,Y ) .	(2.2)
В матричной игре может быть несколько ситуаций равновесия.
Можно доказать следующее свойство.
Пусть	X * Sm , Y * Sn , v – действительное число, тогда для то-
го чтобы	X * ,Y * были оптимальными стратегиями, а v – ценой иг-

ры, необходимо и достаточно, чтобы выполнялись соотношения

100

		n			≤ i ≤ m;
	H (i,Y ) = ∑aij y* j ≤ v, 1
		j=1			(2.3)
		m
			*	i ≥ v, 1	≤ j ≤ n.
	H (X , j) = ∑aij x
		i=1

Под решением матричной игры будем понимать нахождение векторов X и Y , а также значения цены игры v .

2.1.2. Основы метода фон Неймана

Данный метод представляет собой итеративный метод приближенного решения матричных игр размерности ( m ×n ), обладающий достаточно хорошей скоростью сходимости.

Пусть дана платежная матрица

a	a		.	a
11	12			1n
A = .		. .		.	.
am1	am1 .			amn

Положим, что решения в чистых стратегиях нет. Будем искать решение игры (m ×n) в смешанных стратегиях в виде

X = (x , ..., x

)т ,

x

≥ 0,

i =

,

m

i

1, m

1

(2.4)

Y = ( y , ..., y

)т ,

y

≥ 0,

j =

m

j

1, n

,

1

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

m

n

∑xi

=1,

∑y j =1.

i=1

j=1

Тогда, если X * ,Y * – оптимальные смешанные стратегии первой

и второй стороны, то в соответствии с соотношениями (2.3) выполняются неравенства

	n
		≤ v,	i =1, m;
	∑aij y*j
	j=1		(2.5)
	m
	∑aij xi* ≥ v,		j =			,
					1, n
	i=1

где v – цена игры.

101

Суть метода фон Неймана состоит в численном решении системы линейных неравенств (2.5) путем сведения этой задачи к задаче минимизации функции:

			n		m
f (X ,Y ,v) = ∑[z j (X ,v)]2 +∑[ui (Y ,v)]2 ,						(2.6)
			j=1		i=1
где
			при		m
		0		v − ∑aij xi ≤ 0;
z j =					i=1	(2.7)
			m		m
		v − ∑aij xi			при v − ∑aij xi ≥ 0;
			i=1		i=1
			при	n
	0			∑aij x j − v ≤ 0;
ui				i=1		(2.8)
	=	n			n
	∑aij y j − v				при ∑aij y j − vl ≥ 0.
	i=1				i=1

Функция (2.6) представляет собой сумму квадратов невязок правых частей неравенств (2.5).

Точное решение матричной игры соответствует таким значениям аргументов, при которых достигается минимум функции f ( X ,Y , v) , т.е.

f (X * ,Y * , ν) = min f (X ,Y , ν) = 0 .

Задача минимизации функции (2.6) относится к области нелинейного программирования и представляет собой задачу минимизации функции многих переменных с ограничениями. Ограничения наложены на компоненты векторов X и Y в соответствии с определением смешанной стратегии.

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

102

2.1.3. Алгоритм фон Неймана

Алгоритм заключается в следующем.

1. Начальные оценки векторов X и Y и цены игры ν задаются в виде

	xi(0) =	1		,	i =				;
						1, m
		m
	y(j0) =	1		,	j =			;		(2.9)
							1, n
		n

v(0) = max(aij ) + min(aij ) . 2

2. На (l +1) -й итерации новые значения оценок определяются на основании соотношений

l+1

l

l+1

l

~l+1

− θ

), i =1, m ;

xi

= xi

(xi − xi

l+1

l

l+1

l

~l+1

),

j =1, n;

(2.10)

y j

= y j − θ

( y j

− y j

v

l+1

= v

l

− θ

l+1

(v

l

~l+1

),

− v

где

θl+1 =

f (X l ,Y l , vl )

.

f (X

l

,Y

l

, v

l

~ l+1

~l+1

~l+1

)

) + f (X

,Y

, v

~l+1

~l+1

(i =1, m;

j =1, n)

представляют собой оценку

Значения xi

, y j

величин xi* , y*j

на шаге ( l +1). Эти значения рассчитываются с

учетом соотношений (2.7) и (2.8) по формулам:

ul

m

~ l+1

i

при sl

≠ 0,

sl

= ∑uil ;

=

l

(2.11)

xi

s

i=1

l

при

s

l

= 0;

xi

103

(l)

1

(l)

1

(l)

max aij +min aij

tl = 0

~l +1

l

xi

=m

, y j

=n,v

=

2

=

z j

y j

tl

l

n

l

l

∑

−v

ui

=

aij y j

j

=1

~l +1

l

l

l

m

l

=

= v

y j

y j

z j

− ∑ aij xi

i=1

~l+1

n

m

~l

+1~l+1

= ∑

∑ a

l

l

v

x

y

j

j

=1i

=1 ij i

ui

= max(0, ui )

~l+1

n

~l+1

~l+1

l

l

∑

z j

= max(0, z j )

ui

=

j

=1 aij y j

−v

l

m

l

2

n

l

2

~l+1

~l+1

m

~l+1

f

=

∑

(u

)

+

∑

( z

)

−

i

j =1

j

z

j

= v

i

∑ a x

i=1

=1 ij i

~l +1

m

~l +1

)

2

n

~l

+1

)

2

f

= ∑ (ui

+ ∑

( z j

sl

m

f l

ε

i=1

j

=1

= ∑ ul

l

i=1

i

θl

+ 1 =

f

l

~l

+ 1

f

+

f

~l +1

l

s

l

= 0

=

ui

l +1

l

l +1

l

~l +1

x

sl

=

+ θ

)

i

xi

xi

( xi

− xi

l +1

=

l

+ θ

l +1

l

−

~l +1

)

y j

y j

( y j

y j

n

v

l +1

= v

l

+ θ

l +1

(v

l

~l +1

)

~l +1

l

l

l

− v

=

t

=

∑

z

xi

xi

j

j =1

l ≤ N + 1

Рис. 2.1. Структурная схема алгоритма фон Неймана

104

~ l+1 y j

zlj

=t lylj

	n
при t l	≠ 0, t l = ∑zlj ;	(2.12)
	j=1
при t l	= 0.

Оценки цены игры на ( l +1)-шаге вычисляются по формуле

~l+1	~ l+1	~ тl+1	,
v	= X	AY

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

~l+1	m n	~l+1	~l	+1	.	(2.13)
v	= ∑ ∑aij xi		y j

i=1 j=1

3. Итерационный процесс заканчивается, когда значение функции (2.6) становится достаточно малым, т.е. когда выполняется соотношение

f l ≤ ε ,

(2.14)

где ε – заданная точность вычисления.

Структурная схема алгоритма фон Неймана приведена на рис. 2.1.

2.1.4.Математическое программирование

втеории биматричных игр

2.1.4.1.Биматричные игры. Основные теоретические сведения

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

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

Биматричная игра определяется следующим образом. Два игрока обозначены буквами A и B . Игрок A имеет в своем распоряжении m чистых стратегий, а B – n чистых стратегий. Если A выбирает свою i -ю чистую стратегию, а B – j -ю чистую страте-

гию, то выигрыш игрока A есть aij , а выигрыш игрока B есть bij .

105

Определим матрицы A1 и B1 как матрицы размером m × n , (i, j) элементы которых есть aij и bij . Игра полностью определена, ко-

гда заданы матрицы выигрышей A1 и B1 .

Смешанная стратегия игрока A есть столбец X из неотрицательных элементов xi , которые представляют собой относитель-

ную частоту, с которой A будет выбирать свою i -ю чистую стратегию. Таким образом,

x1 + x2 + ... + xm =1.

(2.15)

Аналогично смешанная стратегия игрока B есть столбец Y , неотрицательные элементы y j которого в сумме равны 1. Если A и

B всегда выбирают чистую стратегию случайным образом в соответствии с распределением вероятностей, заданным векторами X и Y , то ожидаемые выигрыши игроков A и B получаются

m n	m n
VA = ∑∑xi aij y j = (X , A1Y ),	VB = ∑∑xibij y j = (X , B1Y ). (2.16)
i=1 j=1	i=1 j=1

Вэтих соотношениях выражение в скобках означает скалярное произведение соответствующих векторов.

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

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

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

( X 0 ,Y 0 ) , такая, что	для любых	других	смешанных	стратегий
(X , Y ) выполняются соотношения
(X 0 , A Y 0 ) ≥ (X , A Y 0 ) , (X 0 , B Y 0 ) ≥ (X 0 , B Y ) .				(2.17)
1	1	1	1

Иными словами, ситуация равновесия – это такая ситуация, отклонение от которой одного из игроков не может увеличить его выигрыш.

106

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

Теорема Нэша гарантирует существование ситуаций равновесия в биматричных играх [40], но не дает никаких средств для их нахождения.

Наиболее эффективным из всех известных алгоритмов для практических вычислений ситуации равновесия является алгоритм, предложенный Лемке и Хоусоном в 1963 г. [59]. По своей сути этот алгоритм тесно связан с методами нелинейного программирования, изложенными в п. 1.3. Следует особо подчеркнуть то обстоятельство, что данный алгоритм может быть без особых затруднений реализован на ЭВМ.

Докажем, прежде всего, две теоремы, на основании которых можно сформулировать постановку задачи.

Теорема 2.1. Ситуация (X 0 ,Y 0 ) является ситуацией равновесия в биматричной игре с матрицами выигрышей A1 и B1 в том и только в том случае, когда

(X 0 , A Y 0 )L	m	≥ A Y 0 ,
1		1		(2.18)
(X 0 , B Y 0 )L		≥ BT	X 0
	n			,
1		1

где Lm и Ln – векторы размерности m и n соответственно, со-

ставленные из единиц (знаки неравенств используются для покоординатного сравнения векторов).

Доказательство. Пусть (X 0 ,Y 0 ) – ситуация равновесия, т.е.

(X 0 , A Y 0 ) ≥ (X , A Y 0 ) при всех X ≥ 0 , (X , L				m	) =1	,
1	1
(X 0 , B Y 0 ) ≥ (X 0	, B Y ) при всех Y ≥ 0 , (Y , L	n	) =1.			(2.19)
1	1

Возьмем в качестве X и Y векторы, у которых одна из компонент равна единице, а остальные – нулю. Получим m + n соотношений

(X 0 , A Y 0 )L	m	≥ A Y 0 ,
1		1		(2.20)
(X 0 , B Y 0 )L		≥ B т	X 0
	n			,
1		1

Наоборот, пусть выполняются указанные условия.

107

единиц. Вследствие процедуры Теорема 2.2. Ситуация

Вектор

X = (x1 , x2 ,..., xm ) , (X , Lm ) =1

(2.21)

является произвольным вектором. Представим его в базисе единичных векторов f1, f2 ,..., fm следующим образом:

X = x1 f1 + x2 f2 +... + xm fm .										(2.22)
Тогда
(X , A Y 0 ) = x ( f , A Y 0 ) +... + x						m	( f	m	, A Y 0 ) ≤
1		1	1 1						1
≤ x ( X 0	, A Y 0 ) + ... + x			m	( X	0 , A Y 0 ) =
1		1						1
= (x + ... + x		m	) ( X 0 , A Y 0 ) = ( X 0 , A Y 0 ) ,							(2.23)
1			1						1

что и требовалось доказать.

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

да к платежным матрицам A и	B , элементы которых – положи-
тельные числа.
Рассмотрим матрицы
A = dE − A1 ,	B = dE − B1 ,	(2.24)

где d – достаточно большая положительная константа, такая, что A > 0 и B > 0 , E – матрица размерности (m × n) , составленная из

(2.24) изменяются и условия (2.20).

X

*

X 0

=

,

(X

* , Lm )

(2.25)

Y *

Y

0

=

(Y

*

, Ln )

является ситуацией равновесия в игре с матрицами выигрышей A1

и

B в том и только в том случае, если пара (X * ,Y * ) удовлетворя-

1

ет соотношениям

Bт X * ≥ Ln , X * ≥ 0 ,

(Y * , Bт X * − Ln ) = 0 ,

(2.26)

108

			AY * ≥ Lm ,				Y * ≥ 0 ,						(X * , AY * − Lm ) = 0 .												(2.27)
Доказательство. Сначала покажем, что
				VA = (X 0 , AY 0 ) = d −													1			.					(2.28)
																(Y	* , Ln )
Имеем				1
										(X * , ((dE − A )Y								* − L			m	)) = 0 ,			(2.29)
	(X * , Lm ) (Y * , Ln )														1
т.е.															1
		0		0			0				0									0
(X			, dEY		) − X			,		A Y			−					(X			, L		m	) = 0 .	(2.30)
														(Y * , Ln )
										1
Учитывая, что EY 0 = Lm , а (X 0 , Lm ) =1 , получим
					(X	0 , A Y 0 ) = d −									1				.						(2.31)
								1							(Y * , Ln )

Это соотношение можно использовать при вычислении цены

игры для участника с платежной матрицей A1 .
Поскольку Y * ≥ 0 , т.е. (Y * , Ln ) ≥ 0 , неравенство
	(dE − A )Y * ≥ L							n		(2.32)
можно преобразовать к виду				1
							Lm
	(dE − A )Y 0				≥					(2.33)
				1		(Y * , Ln )
или
				1
d −					L	m	≥	A Y 0 .		(2.34)
			*
		(Y		, Ln )				1

Используя полученный в начале доказательства результат, при-

дем к неравенству
(X 0 , A Y 0 )L	m	≥ A Y 0 .	(2.35)
1		1
По аналогии доказывается справедливость соотношения
(X 0 , B Y 0 )L	n	≥ B т X 0 .	(2.36)
1		1

109

Согласно теореме 2.1 выполнение этих условий является необ-

ходимым и достаточным, чтобы утверждать, что пара			(X 0 ,Y 0 ) –
ситуация равновесия.
Теорема доказана.
Соотношения
( fi , X * ) ((ai ,Y * ) −1) = 0	для	i =1,..., m ,	(2.37)
(e j ,Y * ) ((b j , X * ) −1) = 0	для	j =1,..., n .	(2.38)
где вектор ai – i -я строка матрицы	A , а вектор b j –		j -й столбец

матрицы B , есть эквивалентная форма записи условий теоремы.

2.1.4.2.Нахождение ситуации равновесия

вбиматричных играх

Доказав теоремы 2.1 и 2.2, получили критерий, с помощью которого удобно определять, является ли пара (X 0 ,Y 0 ) ситуацией

равновесия.

Действительно, если исходить только из определения ситуации равновесия, то для проверки пары (X 0 ,Y 0 ) необходимо переби-

рать все векторы X и Y из множества смешанных стратегий. Сделать это невозможно. Вот почему потребовалось найти необходимое и достаточное условие, в записи которого фигурируют только векторы X 0 и Y 0 .

Постановка задачи нахождения ситуации равновесия. Пусть

A и B > 0 . Задача состоит в том, чтобы найти пару (X * ,Y * ) , для которой выполняются условия:

X * ≥ 0 , Y * ≥ 0 .	(2.39)
( fi , X * ) ((ai ,Y * ) −1) = 0 , (ai ,Y * ) ≥1 для i =1,..., m .	(2.40)
(e j ,Y * )((b j , X * ) −1) = 0 , (b j , X * ) ≥1 для j =1,..., n .	(2.41)

Будем называть любую пару (X * ,Y * ) ситуацией равновесия.

Допустимое множество для поиска смешанных стратегий игрока 1 определяется в соответствии с соотношениями (2.40) и (2.41) следующим образом:

110

Χ = {X | X ≥ 0, B т X − Ln ≥ 0},			(2.42)
т.е. многогранник Χ образуется из векторов X , для которых
( fi , X ) ≥ 0	для	i =1,..., m ,	(2.43)
(b j , X ) ≥1	для	j =1,..., n .	(2.44)
Конец вектора X лежит на границе множества Χ , если хотя бы
одно из этих соотношений является равенством.
Аналогично можно рассмотреть множество
Υ = {Y | Y ≥ 0, AY − Lm ≥ 0}.			(2.45)

Граница множества Υ состоит из точек,				удовлетворяющих хотя
бы одному из m + n уравнений
(ai ,Y ) =1		для	i =1,..., m ,		(2.46)
(e j ,Y ) = 0		для	j =1,..., n ,		(2.47)
а для внутренних точек выполняются неравенства
(ai ,Y ) >1		для	i =1,..., m ,		(2.48)
(e j ,Y ) > 0		для	j =1,..., n .		(2.49)
Легко убедиться, что Χ и		Υ – это многогранники, у которых
имеются соответственно m и n бесконечных ребер.
	*	*	*
Пример 2.1. Рассмотрим игру с матрицами выигрышей
1.75	1.5		1.5		1.8
A1 =	1.8	,	B1 =		.
1.66			1.67		1
Выберем d = 2 , получим матрицы
0.25	0.5		0.5	0.2
A =	0.2	,	B =		.
0.34			0.33		1

Множество Υ , показанное на рис. 2.2, определяется следующими соотношениями:

0.25y1 + 0.5y2 −1 ≥ 0,

0.34 y1 + 0.2 y2 −1 ≥ 0,

y1 ≥ 0, y2 ≥ 0,

111

а множество Χ , показанное на рис. 2.3, – соотношениями:
0.5x1 + 0.33x2 −1 ≥ 0,
				−1 ≥ 0,
0.2x1 + y2
x		≥ 0, x	2	≥ 0.
	1
Y2
6
4
					Yвн
Yгр
2
0	2			4	6	Y1

Рис. 2.2. Множество Υ : Yвн – внутренняя точка; Yгр – точка на границе

Заметим, что здесь и в дальнейшем используется важное предположение, сделанное при постановке задачи:

A > 0 и B > 0 .

Именно это условие гарантирует существование выпуклых открытых многогранников Χ и Υ .

112

Рис. 2.3. Множество Χ :

X 0 , X 1 , X 2

– крайние точки; l, r

– открытые ребра

* * *

Пусть для любого X через p(X ) обозначено множество тех из векторов fi или b j , для которых

( fi , X ) = 0 или (b j , X ) =1.

(2.50)

Множество p(X ) может быть единственным образом записано в виде матрицы

p(X ) = ( p1 , p2 ,..., pr ) ,

(2.51)

где pl (l =1,..., r) – это соответствующий определению вектор fi или b j .

113

Предположим, что матрица B удовлетворяет условию невырожденности: пусть столбцы матрицы C являются столбцами из матрицы (I, B) , где I – единичная матрица, тогда если C = p(X )

для некоторого X Χ, то ранг матрицы C равен числу ее столбцов. Рассмотрим следствия, вытекающие из предположения о невы-

рожденности.

Следствие 1. Если для некоторого X Χ множество p(X ) со-

держит	m векторов, то X полностью определяется своими векто-
рами fi	и b j . В этом случае	X	называется крайней точкой мно-
жества Χ .
Следствие 2. Пусть для X 0 Χ матрица p(X 0 )								имеет ранг r ,
матрица	p(X 0 ) есть матрица размером					(m ×r) ( r		может быть и
равным нулю):
	p(X 0 ) =	( p	, p	2	,..., p	r	) .	(2.52)
		1

Ввиду того, что p1, p2 ,..., pr – линейно независимые векторы, можно присоединить к ним векторы pr+1, pr+2 ,..., pm так, чтобы матрица

D = ( p1 , p2 ,..., pm )							(2.53)
была невырожденной.
Существует обратная матрица
(D−1 )т = ( p1, p2 ,..., pm ) .							(2.54)
Это означает, что
( pi , p j ) = δij =	1	при	i = j,				(2.55)
		при	i ≠ j.
	0
Теорема 2.3. Пусть X 0 Χ и p(X 0 ) = ( p , p				2	,..., p	r	) . Сущест-
вует такая постоянная K , что при			1
m
∑λ2i	≤ K ,						(2.56)

i=1

где λi ≥ 0 для 1 ≤ i ≤ r , точка X , определяемая соотношением

114

		m
X = X 0		+ ∑λi pi ,	(2.57)
		i=1
принадлежит множеству Χ .
Доказательство. Пусть	p – любой столбец матрицы (I, B). То-
гда
		m
( p, X ) = ( p, X 0 ) + ∑λi ( p, pi ) .			(2.58)
		i=1
Рассмотрим два случая.
1. Если p p(X 0 ) , т.е.	p = pi , 1 ≤ i ≤ r .

На основании следствия 2 имеем ( p, X ) = ( p, X 0 ) + λi . При этом

1,	если p	− столбец из B,	(2.59)
( p, X 0 ) =	если p	− столбец из I.
0,

По условию теоремы λi ≥ 0 , 1 ≤ i ≤ r , следовательно, ( p, X ) ≥

≥( p, X 0 ) и X является подмножеством Χ .

2.Если p p(X 0 ) , т.е. p – любой другой столбец из (I, B). В этом случае

1,	если p	− столбец из B,	(2.60)
( p, X 0 ) >	если p	− столбец из I.
0,

Пусть p – столбец из B , тогда ( p, X 0 ) =1 + α , где α > 0 . Выбрав λi , 1 ≤ i ≤ m , достаточно малыми, такими, чтобы

	m
	∑λi ( p, pi )	≤ α,	(2.61)
	i=1
получим ( p, X ) ≥1.

По аналогии можно рассмотреть столбец матрицы I . Таким образом, при соответствующем выборе λi , 1 ≤ i ≤ m , из соотношения

( p, X 0 ) ≥ 0 следует, что ( p, X ) ≥ 0 , а из соотношения ( p, X 0 ) ≥1 следует, что ( p, X ) ≥1. Следовательно, X Χ.

Теорема доказана.

115

Таким образом, доказана возможность перемещения вдоль границы многогранника Χ или Υ от одной крайней точки к другой. Именно в таком перемещении и проверке в каждой крайней точке условий (2.37) и (2.38) и заключается суть алгоритма Лемке – Хоусона поиска ситуации равновесия.

Обратимся вновь к рис. 2.3. Точка Х0 является крайней точкой множества Χ . Поскольку выполняются условия

0.5x0	+ 0.33x0		=1	и	0.2x0			+ x		0	=1 ,	(2.62)
1		2					1		2
матрица p(X 0 ) невырождена:
		p(X	0 ) = ( p		, p	2	) ,					(2.63)
где				1
		0.5					0.2
p1				а	p2							(2.64)
		=	,			=		1	.
		0.33
Согласно соотношению (2.57) точки
		X = X 0 + λ1 p1 + λ2 p2										(2.65)

при малых положительных λi принадлежат Χ . Пусть λ2 = 0 . Для

X = X 0 + λ1 p1 ,	λ1 > 0 ,	(2.66)
матрица p(X ) получается из p(X 0 )	вычеркиванием только одно-

го столбца, а именно р1. Такие точки образуют открытое ребро многогранника Х с концом Х0. На рис. 2.3 оно обозначено буквой r.

Для X 1 Χ матрица p(X 1 ) = ( p ) .	При достаточно малом λ	2
1
точки
X = X 1 + λ2 p2	(2.67)

принадлежат множеству Χ и для них p(X ) = p(X 1 ) . Такие точки

образуют открытое ребро 1, содержащее X 1 . В общем случае матрица p(X ) , где X принадлежит открытому ребру Χ , имеет ранг,

равный (m −1) .

На рис. 2.3 показаны два неограниченных ребра, каждое из которых обладает лишь одной концевой точкой X = kfi с соответствующим k . Таким образом, у многогранника Χ имеется ровно m

116

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

смежных крайних точек,	например, X 2			и	X 0	отличаются лишь
одним столбцом. В нашем примере
0.5	1		, а p(X 0 )	0.5		0.2		(2.68)
p(X 2 ) =	0			=	0.33	1	.
0.33

Аналогичные рассуждения можно провести для множества Y. Пусть q(Y ) – матрица, подобная рассматривавшейся для эле-

ментов X множества Χ матрице p(X ) . Предположим также, что

матрица A удовлетворяет условию невырожденности, аналогично тому, которое было наложено на B.

Рассмотрим Ζ = Χ×Υ . Точка Z = (X ,Y ) из Ζ называется край-

ней для Z, если X – крайняя точка многогранника X, а Y – крайняя точка Y. Очевидно, что такая точка должна удовлетворять (m + n)

условиям из системы

( fi , X ) = 0,				i =1,..., m,
			= 0,	j =1,..., n,
(b j , X ) −1					(2.69)
				j =1,..., n,
(e j ,Y ) = 0,
(a		,Y ) −1 = 0,		i =1,..., m.
	i

Будем говорить, что Z = (X ,Y ) лежит на открытом ребре много-

гранника Ζ , если одна из ее координат является крайней точкой Χ или Υ , а другая лежит на открытом ребре. Для точек, принадлежащих открытому ребру, выполняются (m + n −1) из указанных

выше соотношений.

Теорема 2.4. Любая ситуация равновесия для невырожденной задачи является крайней точкой множества Ζ .

Доказательство.	Если	Z * = (X * ,Y * )	удовлетворяет условиям
постановки задачи,	то	для i =1,..., m	или	( fi , X * ) = 0 ,	или
(ai ,Y * ) −1 = 0 , и	для	j =1,..., n	или	(e j ,Y * ) = 0 ,	или

(b j , X * ) −1 = 0 . Из условия невырожденности следует, что X Χ может удовлетворять не более чем m соотношениям вида:

117

( fi , X * ) = 0 или (b j , X * ) −1 = 0 ,

(2.70)

а Y Υ может удовлетворять не более чем n соотношениям вида:

(e j ,Y * ) = 0 или (ai ,Y * ) −1 = 0 .	(2.71)
Но (X * ,Y * ) должны удовлетворять, по меньшей мере,	(m + n)

таким соотношениям. Таким образом, удовлетворяются в точности

(m + n) соотношений, т.е. для X *			выполнены ровно m из соотно-
шений	( fi , X ) = 0,		i =1,..., m,
					(2.72)
	(b j , X ) −1 = 0,		j =1,..., n.
Следовательно,	X * – крайняя точка Χ . Аналогично Y *				– край-
няя точка Υ , а Z * = (X * ,Y * )		– крайняя точка Z, что и требовалось
доказать.
Одновременно	получен	следующий		критерий:	если

Z* = (X * ,Y * ) – ситуация равновесия, то для любого s, 1 ≤ s ≤ m + n , или s -й столбец матрицы (I, B) принадлежит p(X * ) , или s-й

столбец матрицы (Aт , I ) принадлежит q(Y * ) , но не тот и другой

одновременно. Подчеркнем, что это утверждение имеет силу только для ситуаций равновесия.

2.1.4.3. Метод Лемке – Хоусона. Теоретические основы

Механизм алгоритма Лемке – Хоусона заключается в том, чтобы, последовательно двигаясь от одной крайней точки множества Ζ к другой крайней точке, за конечное число шагов найти одну из ситуаций равновесия. Необходимо доказать несколько теорем, на основании которых строится схема движения. Прежде всего нужно получить само множество путей движения. Для этого рассмотрим точки Z = (X ,Y ) Ζ , удовлетворяющие хотя бы (m + n −1) из ус-

ловий

( fi , X * ) ((ai ,Y * ) −1) = 0 ,	i =1,..., m ,
(e j ,Y * ) ((b j , X * ) −1) = 0 ,	j =1,..., n .	(2.73)

118

Пусть через H s обозначено множество точек		Z Ζ , удовле-
творяющих всем этим уравнениям, кроме, возможно, уравнения
(es ,Y * ) ((bs , X * ) −1) = 0 .		(2.74)
Пример 2.2. Построить H1 для биматричной игры, рассмотрен-
ной в примере 2.1. Для Z = (X ,Y ) H1	одновременно должны вы-
полняться соотношения:
y2 = 0 или 0.2x1 + x2 −1 = 0	(прямая 4),
x1 = 0 или 0.25y1 + 0.5y2 −1 = 0	(прямая 1),
x2 = 0 или 0.34 y1 + 0.2 y2 −1 = 0	(прямая 2),
а соотношение
y1 = 0 или 0.5x1 + 0.33x2 −1 = 0	(прямая 3)

может не выполняться.

На рис. 2.4 и 2.5 показаны точки (X ,Y ) (X Χ,Y Υ), принадлежащие множеству H1 . Множество H1 состоит из:

– открытого ребра с концом в крайней точке (X 0 ,Y 0 ) , образованного точками (X ,Y 0 ) , где X удовлетворяет условию

0.2x1 + x2 −1 = 0 ;

–открытого ребра с концом в крайней точке (X 1,Y 1 ) , образо-

ванного точками (X 1 ,Y ) , где Y удовлетворяет условию

0.25y1 + 0.5y2 −1 = 0 ;

– неограниченного ребра (X ,Y 1 ) с концом (X 1,Y 1 ) .

Полученный результат можно обобщить, доказав теоремы 2.5 и 2.6.

Теорема 2.5. Любая точка из H s или является крайней для Ζ ,

или лежит на открытом ребре Ζ .

Доказательство. Если Z H s удовлетворяет всем (m + n) со-

отношениям, то она является крайней. Если же она удовлетворяет
(m + n −1) из уравнений
( fi , X ) ((ai ,Y ) −1) = 0 ,	i =1,..., m ,
(e j ,Y ) ((b j , X ) −1) = 0 ,	j =1,..., n ,	(2.75)

119

	Рис. 2.4. Проекция s − пути на плоскость Y
X2
6
		Многогранник X
4	X1=0
	X2
2	3
	X0	4
0		X1 X2=0		X2
	2	4	6	X1
	Рис. 2.5. Проекция s − пути на плоскость Χ
		120

то в силу условий невырожденности выполнены также (m + n −1) из соотношений

( fi , X ) = 0,

i =1,..., m,

j =1,..., n,

(b j , X ) −1 = 0,

(2.76)

j =1,..., n,

(e j ,Y ) = 0,

(a

,Y ) −1 = 0,

i =1,..., m.

i

Следовательно, эта точка лежит на ребре.

Теорема 2.6. Точка множества H s

образует единственное неог-

раниченное ребро многогранника Ζ .

Доказательство. Рассмотрим множество

H

1

. Пусть Y 1 = k

0

e .

1

При соответствующем k0 точка

Y 1

будет крайней для Υ (см.

рис. 2.4). Так как (e j ,Y 1 ) = 0 при

j ≠1 и этих соотношений всегда

(n −1), то так как Y 1

– крайняя точка Υ , следует, что она должна

быть

определена

еще

одним

соотношением,

например

(ar ,Y 1 ) −1 = 0 . Итак,

q(Y 1 ) = (e2 ,..., en , ar ) .

(2.77)

Для случая, рассмотренного в примере 2,

r =1 . При достаточно

больших

k точки X = kfr

принадлежат

Χ .

Пусть k1

таково, что

X 1 = k

f

r

является крайней точкой Χ . Поскольку соотношения

1

(e j ,Y 1 ) ((b j , X ) −1) = 0 ,

j = 2,..., n ,

( fi , X ) ((ai ,Y 1 ) −1) = 0 ,

i =1,..., m ,

(2.78)

где X = kfr , (k ≥ k1 ) , выполняются, то можно утверждать, что пары (X ,Y 1 ) составляют открытое ребро многогранника Ζ , лежащее

в H1 . Точка (X 1 ,Y 1 ) является концом этого ребра.

Докажем единственность этого ребра. Рассмотрим любое другое неограниченное ребро многогранника Ζ , например (X ,Y 0 ) , где X = kfl , (l ≠ r) , а Y 0 – крайняя точка множества Υ .

121

Такое ребро не принадлежит H1 в силу того, что либо условие

( fl , X ) ((al ,Y 0 ) −1) = 0 ,	(2.79)
либо одно из условий
(e j ,Y 0 ) ((b j , X ) −1) = 0 , j = 2,..., n ,	(2.80)

не выполняется. Теорема доказана.

Два открытых ребра многогранника Ζ называются смежными, если они имеют общий конец. Последовательность смежных открытых ребер из H s вместе с их концевыми точками называется

s -путем.

При решении задачи, рассмотренной в примере 2.2 (см. рис. 2.4 и 2.5), можно двигаться по следующему s -пути с началом в точке

(X 2 ,Y 1 ) :

– двигаемся вдоль неограниченного ребра (X ,Y 1 ) до точки

(X 1,Y 1 ) ;

–из точки (X 1,Y 1 ) продолжаем движение вдоль ребра (X 1 ,Y ) , то есть по прямой 1 в множестве Υ , попадаем в крайнюю точку

(X 1,Y 0 ) ;

–(X 1,Y ) является концом ребра (X ,Y 0 ) , где точки X Χ

принадлежат прямой 4, по этому ребру можно попасть в (X 0 ,Y 0 )

– конечную точку рассматриваемого s -пути.

Теорема 2.7. Пусть Z – крайняя точка Ζ и Z H s . Найдутся одно или два открытых ребра, лежащих в H s и имеющих Z своей

крайней точкой; Z является ситуацией равновесия в том и только в том случае, когда такое ребро единственное.

Доказательство. Положим, что Z – крайняя точка Ζ и Z H s . Возможны два случая.

1. Пусть (es ,Y ) ((bs , X ) −1) = 0 . Поскольку справедливы (m + n) соотношений

( fi , X ) ((ai ,Y ) −1) = 0 , i =1,..., m ,

122

(e j ,Y ) ((b j , X ) −1) = 0 , j =1,..., n ,

(2.81)

точка Z = (X , Y ) – ситуация равновесия. Как было показано, при этом или (es ,Y ) = 0 , или (bs , X ) −1 = 0 , но не оба сомножителя равны нулю сразу.

	Допустим, что (bs , X ) −1 = 0		( (es ,Y ) = 0		рассматривается ана-
	логично). Так как Z = (X ,Y )	– крайняя точка Ζ , выполнено m из
	соотношений
	( fi , X ) = 0,			i =1,..., m,	(2.82)
	(b j , X ) −1 =		0,	j =1,..., n,
	и n из соотношений
	(e j ,Y ) = 0,			j =1,..., n,	(2.83)
	(ai ,Y ) −1 =		0,	i =1,..., m.
	Следовательно, существует	(m + n) ребер			многогранника Ζ ,

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

нарушается условие (bs , X ) −1 = 0 .

2. Пусть (es ,Y ) ((bs , X ) −1) > 0 . Так как Z = (X , Y ) –	крайняя
точка Ζ , X определяется m векторами из множества
{f1,..., fm ,b1,...bn },	(2.84)
а Y определяется n векторами множества
{e1,..., en , a1,..., am }.	(2.85)
В нашем случае (es ,Y ) > 0 и (bs , X ) −1 > 0 . Но Z H s ,	поэтому

найдется такой индекс q , для которого выполнены сразу или оба соотношения

( fq , X ) = 0	и	(aq ,Y ) −1 = 0 ,		(2.86)
или оба соотношения
(eq ,Y ) = 0	и	(bq , X ) −1 = 0 .		(2.87)
Рассмотрим первый вариант.		В	H s лежат два ребра –	одно,
вдоль которого нарушается условие			( fq , X ) = 0 , и другое,	вдоль

которого нарушается условие (aq ,Y ) −1 = 0 . При этом Z не явля-

123

ется ситуацией равновесия. На любом другом ребре Ζ , имеющем Z концевой точкой, не выполнено какое - либо из соотношений, определяющих H s .

Теорема доказана.

Доказанные выше теоремы гарантируют существование крайней точки Z , лежащей в H s . Начав с такой точки Z1 , можно двигаться

по ребру, лежащему в H s : или это ребро закончится в другой крайней точке Z2 , или это единственное неограниченное ребро, лежащее в H s . В первом случае или Z2 – ситуация равновесия, и процесс не может быть продолжен, или существует другое ребро с концом в Z2 , лежащее в H s , по которому можно продолжать двигаться. Процесс прекращается в следующих случаях, когда:

–попадаем на неограниченное ребро;

–достигаем ситуации равновесия, отличной от точки Z1 ;

–приходим в точку, в которой уже были.

Вернуться можно только в начальную точку пути, так как в противном случае имелась бы точка, к которой примыкают три ребра. Таким образом, начав с Z1 , возвращаемся в Z1 , или нет. В первом

случае путь называется замкнутым (рис. 2.6, а). Если путь не замкнут, то он заканчивается либо в ситуации равновесия (рис. 2.6, б), либо на неограниченном ребре (рис. 2.6, в).

Рис. 2.6. Возможные s-пути

124

Теорема 2.8. Пусть P – s -путь, содержащий неограниченное ребро F . Тогда на P лежит в точности одна ситуация равновесия. Ее можно вычислить, последовательно проходя путь P , начиная с ребра F . Общее число ситуаций равновесия конечно и нечетно.

Доказательство. Единственный s -путь P , начинающийся с F , должен закончиться в некоторой крайней точке, являющейся ситуацией равновесия (см. рис. 2.6, в). Отличный от P незамкнутый путь должен иметь две конечные точки (см. рис. 2.6, б), каждая из которых является ситуацией равновесия. Отсюда следует утверждение теоремы.

2.1.4.4. Алгоритм Лемке – Хоусона решения биматричных игр

Дана биматричная игра с матрицами

A1 и

A2 размерностью

m ×n .

N =

1. Положим

0

– номер

итерации.

Выбираем

d = max (

aij1

,

aij2

)+1.

Определяем матрицы

A = dE − A1

и

i, j

B = dE − A2 .

2. Выбираем

начальные

базисы

q0 (Y ) ={e ,..., e

n

}

и

1

p0 (X ) ={ f1 ,..., fm } . Составляем симплекс-таблицы.

Таблица A

Базис		a1	a2	...	am		e1	e2	...	en
e1		a11	a21	...	am1		1	0	...	0
e2		a12	a22	...	am2		0	1	...	0
...		...	... ... ...				... ... ... ...
en		a1n	a2n	...	amn		0	0	...	1
												Таблица B
Базис		b1	b2	...	bn		f1	f2	...	fm
f1		b11	b12	...	b1n		1	0	...	0
f2		b21	b22	...	b2n		0	1	...	0
...		...	... ... ...				... ... ... ...
fm		bm1	bm2	...	bmn		0	0	...	1

125

3. Выбираем множество H s , в котором будет производиться поиск ситуации равновесия. Определяем исходную точку

(X 0 ,Y 0 ) :

Y 0

=

0, 0, ..., 0,

1 l , 0, ..., 0

,

s-япозиция

где l = a *

s

= min(a1s ,a2s ,...,ams ) , т.е. минимум в s -й строке табли-

i

i

цы A;

X 0

0, 0, ..., 0,

1 k,

,

=

0, ..., 0

i*-япозиция

где k = b

*

j

* = min

(b

*

,b *

2

,...,b *

), т.е. минимум в строке с номе-

i

j

i

1

i

i n

ром i* таблицы B.

Таким образом,

(X 0 ,Y 0 ) выбрана так,

чтобы она была крайней

точкой множества Ζ . Согласно теореме 2.6 (X 0 ,Y 0 ) H s , где H s состоит из точек множества Χ ×Υ , удовлетворяющих уравнениям

( fi , X )((ai ,Y ) −1) = 0 , i =1,..., m ,

(e j ,Y )((b j , X ) −1) = 0 , j =1,..., n , j ≠ s ,

и, возможно, уравнению (es ,Y * )((bs , X * ) −1) = 0 .

4. Если последнее уравнение выполнено для точки (X 0 ,Y 0 ) , то

это – ситуация равновесия. Переходим к шагу 13. В противном случае переходим к следующему шагу.

5. Изменяем базисы p0 и q0 , для этого найдем множества p(X 0 ) ={ fi , b j | ( fi , X 0 ) = 0, (b j , X 0 ) =1} ,

q(Y 0 ) ={e j , ai | (e j ,Y 0 ) = 0, (ai ,Y 0 ) =1}.

126

6. Выписываем таблицы AN и B N ( N – номер итерации) для новых базисов так же, как в симплекс-методе []. Чтобы получить

таблицу AN , надо исключить из базиса вектор es и ввести вектор ai* . Для этого из всех строк таблицы необходимо вычесть s-ю

строку, умноженную на множители, подобранные таким образом, чтобы во всех строках, кроме s-й, i*-й элемент обратился в нуль. Затем разделить s-ю строку на ai*s .

Таблица AN

	Базис	a1	a2	...	am	e1	e2	...	en	λ
		α11	α21	...	αm1	q11	q12	...	q1n	λ1
	q(Y N )	α12	α22	...	αm2	q21	q22	...	q2n	λ2
		...	... ... ...			... ... ... ...				...
		α1n	α2n	...	αmn	qn1	qn2	...	qnn	λn
		ξ1 −1	ξ2 −1	...	ξm −1	y N	y N	...	y N
						1	2		n

В таблицу AN добавляется снизу строка и справа один столбец. Значения ξi = (ai ,Y 0 ) , причем ai берется из таблицы, полученной на предыдущем шаге вычислений; y0j – j -я компонента вектора

Y 0 ; заметим, что

ξi

≥1 , i =1,..., m , так как Y 0 Υ . Значение λ j

для j -й строки получается следующим образом:

1) вычисляются

*

ξk −1

N

yr

λ

=

min

−

,−

≥ 0 ,

α

q

αkj <0,q jr <0

kj

jr

k =1,m,r=1,n

**

ξk −1

N

yr

λ =

min

−

,−

≤ 0

;

α

q

αkj >0,q jr >0

kj

jr

k =1,m,r=1,n

2) так как (X 0 ,Y 0 ) – крайняя точка, можно показать, что либо λ* , либо λ** равно нулю. В столбец λ в качестве λ j добавляет-

127

ся то из чисел λ* и λ** , которое отлично от нуля. Если они оба равны нулю, добавляется нуль.

Аналогично формируется таблица B N . Из базиса исключается fi* и на его место вводится b j* . В нижней строке таблицы поме-

щаются числа η j = (b j , X N )

без единицы и компоненты xiN . Чис-

ла μi рассчитываются аналогично λ j .

Таблица B N

Базис

b1

b2 ...

bn

f1

f2

...

fm

μ

β11

β12 ...

β1n

p11

p12

...

p1m

μ1

p(X N )

β21

β22 ...

β2n

p21

p22

...

p2m

μ2

... ... ... ...

... ... ... ...

...

βm1

βm2 ...

βmn

pm1

pm2

...

pmm

μm

η1 −1

η2 −1 ...

ηn −1

x N

x N

...

x N

1

2

m

7. Утверждаем, что Y = Y 0 + λ j q j , где

q j

– j -я строка входя-

щей в таблицу матрицы

qij

, является крайней точкой множества

Υ . Покажем, прежде всего, что Y Υ , т.е. что

(ak ,Y 0 + λ j q j ) ≥1 .

Матрица

qij

является обратной матрице, столбцами которой слу-

жат векторы базиса, поэтому справедливы соотношения

(ak , q j ) = αkj ,

k =1,2,..., m ,

j =1,2,..., n .

Таким образом,

(ak ,Y 0 ) + λk (ak , q j ) = ξk + λ j αkj ,

причем правую часть равенства преобразуем следующим образом:

ξ

k

−1

(ξ

k

−1) + λ

j

α

kj

= −α

−

− λ

.

α

kj

j

kj

128

	Если λ j > 0 , то −		ξk	−1			*
						− λ j ≥ 0 по определению для	λ , и пото-
			αkj
	му при αkj	< 0 имеем
				(ξk −1)+ λ j αkj ≥ 0 ,
	т.е.
				(ak ,Y 0 ) + λk (ak , q j ) ≥1.
	Если же λ j	> 0 и αkj		> 0 , то λ j αkj > 0 , поэтому вновь
				(ak ,Y 0 ) + λk (ak , q j ) ≥1,
	т.е. доказали, что Y Υ .
	При λ j	< 0 – аргументация аналогична. Таким образом, Y Υ .

Точка Y определена так, что она является крайней точкой множества Υ . Ввиду того, что Y – конец ребра, другим концом которого

является Y 0 , числа λ* и λ** не могут одновременно отличаться от нуля, так как в противном случае

(Y 0 + λ*q j ) Υ , (Y 0 + λ**q j ) Υ , λ* > 0 , λ** < 0 ,

вопреки тому, что Y 0 – крайняя точка.

Таким образом, на данном шаге алгоритма рассчитываем край-

нее точки
X i = X 0 + μi pi ,	i =1,2,..., m ,
Y j = Y 0 + λ j q j ,	j =1,2,..., n .

Точки (X 0 ,Y j ) и (X i ,Y 0 ) являются смежными с точкой

(X 0 ,Y 0 ) .

8.Находим q(Y j ) и p(X i ) в соответствии с определением базиса (п. 5 настоящего алгоритма).

9.Для каждой пары (X 0 ,Y j ) и (X i ,Y 0 ) формируем множество

M :

129

	, f		er		M , еслиer		q(Y	j	),br	p(X
M (X i ,Y j ) = e
r		k	f		M , еслиf		p(X i ), a				q(Y
				k		k				k

10. Если для некоторой пары (X i ,Y j )

M (X i ,Y j ) ={e1, e2 ,..., en , f1 ,..., fm } ,

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

(X * ,Y * ) = (X i ,Y j ) .

i), j ).

(X i ,Y j ) –

Переходим к шагу 13.

Если такой пары не найдется, то переходим к следующему шагу. 11. Выбираем точку Z = (X ,Y ) , для которой в M (X ,Y ) недос-

тает вектора es . Тогда (X ,Y ) H s . По теореме 2.7 найдется пара таких точек Z1, Z2 , являющихся концами ребер, которые пересе-

каются в (X 0 ,Y 0 ) . Если к точке (X 0 ,Y 0 ) примыкает неограниченное ребро, такая точка будет одна.

12. Выбираем ту из найденных точек Z1, Z2 , которая не была

включена в s -путь на предыдущих итерациях, изменяем базис p(X ) или q(Y ) и возвращаемся на шаг 6. (При этом на шаге 6 ме-

няется лишь одна из таблиц A или B .)

Согласно теореме 2.8 через конечное число шагов процесс закончится или в ситуации равновесия, или в исходной точке

(X 0 ,Y 0 ) , если произойдет зацикливание. В последнем случае выбираем новое s и начинаем поиск сначала.

13. Нормируем пару (X * ,Y * ) :
X 0 =	X *		,		Y 0 =		Y *	.
	(X * , Lm )
						(Y * , Ln )
Вычисляем выигрыши игроков:
для игрока 1					1
	V1 = d −						,
				(Y * , Ln )

130

для игрока 2

	V2 = d −	1	.
		(X * , Lm )

(Последние соотношения выводятся при доказательстве теоре-

мы 2.2.)

2.2.Оптимальное распределение нагрузки

всистеме ядерных реакторов

2.2.1.Физическая постановка задачи

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

При этом решение проблемы маневренности энергоблоков АЭС возможно только после проведения комплекса научно-исследова- тельских и опытно-конструкторских работ. Основные из них следующие:

1)снятие ограничений на число циклов пуск-остановка для всего оборудования АЭС;

2)улучшение схем пусков и остановок АЭС и снижение потерь тепла при расхолаживании блоков за счет рационального использования остаточного энерговыделения ядерного топлива и теплоаккумулирующей способности графитовой кладки;

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

4)приведение в соответствие с требованиями эксплуатации систем управления и защиты реакторов и обеспечение реакторов необходимым запасом реактивности.

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

131

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

Q = aρ ,

где a – темп выгорания, 1кВт сут; ρ – запас реактивности, отн. ед.; Q = Qm −Q – разность между энерговыработкой при работе в базисном режиме ( Qm ) и режиме переменных нагрузок ( Q ), кВт сут.

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

снижения мощности реактора ε = W , где Wн – номинальная мощ-

Wн

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

В рамках точечной модели запас реактивности, обеспечивающий работу реактора на пониженной мощности εWн в течение произвольного времени, после снижения мощности равен

ρ = xм − xр ,

где xм и xр – максимальная и равновесная концентрации ксенона,

	нормированные на	ν f Σ f	;	Σ f	– макроскопическое сечение деле-
		σX

ния активной зоны реактора; ν f – среднее число вторичных ней-

тронов на акт деления; σX – сечение поглощения ксенона, см2.

На АЭС, как правило, устанавливается несколько энергоблоков. В общем случае энергоблоки могут отличаться электрической мощностью, темпом выгорания, другими характеристиками и рабо-

132

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

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

2.2.2. Математическая постановка и решение оптимизационной задачи

Математически задача формулируется следующим образом: найти

			N	ρ	i
min			∑
				ai
ρ1 ...		ρN i=1
при ограничениях:
N
∑δi εi (	ρi ) = α;
i=1
0 ≤ ρi	≤		ρi m , i =1, ..., N,

где N – число реакторов на станции; Δρi – оперативный запас реактивности i-го реактора; εi – степень снижения мощности i-го реак-

N

∑Wi

тора; α – заданная степень снижения мощности АЭС, α = Ni=1 ;

∑Wi н

i=1

δi – доля электрической мощности i-го реактора, δi = NWi н .

∑Wi н

i=1

Заметим, что оптимизация возможна только при < 0 < α <1 , так как при α = 0 и α =1 значения оперативных запасов реактивности определены:

133

ρi = ρi м,

ρi = 0.

Поэтому качественно понятно, что эффект от оптимизации будет сильнее всего сказываться в средней части диапазона изменения α . Вид зависимости ε( ρ) возможной степени снижения

мощности реактора от запаса реактивности, вообще говоря, определяется режимом изменения мощности реактора.

Наиболее простым режимом для анализа является снижение мощности реактора с максимально возможной скоростью до определенного уровня и поддержание реактора на данном уровне до момента выхода АЭС на номинальную мощность. Характер зависимости ε( ρ) при различных плотностях потоков нейтронов в

этом случае показан на рис. 2.7.

Возможная степень снижения мощности,ε

	ϕ=1. 1013H/см2с
1,0	ϕ=3. 1013H/см2с
0,8	ϕ=5. 1013H/см2с
0,6
0,4
0,2

0,0

0

1

2

3

4

Величина резервируемого запаса пеактивности, Δρ%

Рис. 2.7. Зависимость возможной степени снижения мощности от величины резервируемого запаса реактивности

134