Математическая экономика

Будем исходить из того, что в нашем распоряжении имеется ЦВМ, позволяющая вычислять значения функции f (а возможно, и ее производных) в любой выбранной пробной точке xk . Решение задачи чис-

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

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

ние искомой величины x*.

Для рассматриваемой одномерной задачи этот подход можно реализовать путем разбиения отрезка [a, b] на N достаточно мелких частей длины h; вычисления функции в полученных точках:

x0 = a; x1 = a + h;...; xk = a + kh;...;

xN = a + Nh = b ; сравнения значений функции y0 , y1,... в этих точках и вы-

бора точки, в которой значение этой функции оказалось наименьшим

(рис. 5.3).

Достоинством этого метода, называемого методом равномерного пассивного поиска, является то, что при достаточно мелком шаге h он позволяет находить решение и в том случае, когда функция имеет несколько экстремумов. Его недостаток – большой объем вычислений, необходимых для отыскания решения с заданной точно-

стью ε (нужно выбирать N ≥ b −ε a ).

Блок-схема алгоритма, реализующего этот метод, представлена на рис. 5.4.

Программирование конкретной функции в этой схеме предусмотрено в блоке 4.

y

аb х

х*sup	х*inf хк	h	хN

Рис. 5.3

Начало

1

Ввод А,В,Е

2H:=E

X:=A

I:=1

4

F:=F(X)

		I>1			Да
	5
Нет					Да	3
						F<FM

6 FM:=F

XM:=X Нет

7X:=X+H

I:=I+1

Да 8

X≤B

Нет

9

Вывод

XM,FM

Конец

Рис. 5.4

141

В тех случаях, когда в рассматриваемом промежутке заданная функция имеет один экстремум, она называется унимодальной. Для минимизации такой функции существует ряд эффективных алгоритмов (метод Фибоначчи, методы золотого сечения, дихотомии, квадратичной аппроксимации и др.) [2, 15, 26, 37]. Рассмотрим один из них – простой и весьма эффективный метод золотого сечения (МЗС). Он основывается на очевидном свойстве унимодальной функции (рис. 5.5, а), в соответствии с которым наличие значений функции в двух точках, выбранных на интервале L, содержащем экстремум и называемом интервалом неопределенности, позволяет отбросить часть этого интервала как заведомо не содержащую экстремального значения, и тем самым сократить интервал неопределенности

(рис. 5.5, б, в).

уy

y

y

y1

у22

у22

yу1

xх

у1

x1

x 2

xх11

x22

х

х2

L

х

а)

L

L

L

б)

Рис. 5.5

Эту процедуру нужно повторять до тех пор, пока длина интервала не

станет меньше допустимой погрешности ε.

С

D

Lk

Особенностью МЗС является то,

Lk

что при выборе пробных точек исполь-

А

В

Lkk+1+

1

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

чек в Lk +1 совпадает с одной из точек в

Рис. 5.6	Lk	(рис. 5.6).

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

ADAB = DBAD ,

называемому золотым сечением. Такое сечение дают точки: xC = xA +(1−τ)(xB − xA), xD = xA + τ(xB − xA),

где τ = (−1+ 5) / 2 = 0,618034 .

142

Блок-схема алгоритма, реализующего МЗС, представлена на рис. 5.7. В блоках 4, 10, 13 предусмотрено обращение к подпрограмме F вычисления заданной целевой функции; величина D представляет собой длину текущего интервала неопределенности.

Начало

1

Ввод А,В,Е

2

Конец 3

15

	Вы-
		4
	вод
	X,Y

14
X:=(X1+X2)/2	5	D:=D*R2
Y:=F(X)

	Да	7		Нет
		Y1<Y2

									11
	8
	X2:=X1								X1:=X2
X1:=X2-D*R									X2:=X1+D*R
		9
										12
						Y2:=Y1								Y1:=Y2
			10								13
						Y1:=F(X1)								Y2:=F(X2)

6

D<E

Да

Нет

Рис. 5.7

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

Рафсона. В соответствии с ним для отыскания точки x* минимума f (x)

143

(5.12)

– произвольная на-

строится числовая последовательность			x0 , x1, x2 ,... по рекуррентному со-
отношению
	x	= x	−	f '(xk )	,	(5.10)
	k +1	k		f "(xk )
где x0 – некоторое начальное значение.
Как только для	некоторого	k будет получено xk , для				которого
f '(xk ) – мало, процесс вычислений прекращают, полагая x* ≈ xk .
Если вычисление	f " невозможно или затруднено, то поиск можно
вести по соотношению
	xk +1 = xk − γ f '(xk ),					(5.11)

в котором γ – некоторое число, подлежащее выбору.

§5.3. Численные методы многомерной безусловной оптимизации

сиспользованием производных

Для поиска	значений xi = xi* (i =		) , при которых ЦФ
		1, n
y = f (x1,..., xn ) = f (x)	принимает наименьшее значение, можно использо-

вать распространение метода Ньютона-Рафсона (5.10) на многомерный случай

x (k +1) = x (k ) − Гk f (x (k ) ),

где Гk = Hk−1; Hk = 2 f (x (k) ) – матрица Гессе; x(0) чальная точка.

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

x(0) , x(1) ,... сходится к искомой экстремальной точке, причем в некоторых случаях скорость сходимости оказывается очень высокой. Так, для квадра-

тичной ЦФ решение при произвольной x(0) получается за один шаг [26]. При практическом использовании (5.12) могут возникнуть сложности с вычислением вторых производных, а главное – с обращением матри-

цы H, так как она часто оказывается близкой к вырожденной.

144

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

	γ	0		0
	1	γ2
Г =	0			0		,
		.	.	.
.
	0	0		γn

где γi – некоторые числа. Тогда в развернутом виде (5.12) принимает следующий вид:

x(k +1)

= x(k)

− γ

∂f

;

1

1

1

∂x1

x =x(k )

.

.

.

.

.

(5.13)

x(k +1)

= x(k)

− γ

∂f

.

n ∂x

n

n

(k )

n

x

=x

Характерно, что при

γ

=... = γ

n

в этом случае переход от x(k) к

1

x(k +1) осуществляется в направлении антиградиента, т.е. направлении наибольшего убывания ЦФ. Поэтому рассматриваемую процедуру поиска минимума называют методом градиентного спуска. Возможны различные модификации этого метода, отличающиеся способом выбора длины шага в указанном направлении. Так, при градиентном методе с дроблением шага [37] числа γ1,..., γk принимаются равными некоторому числу γ = γ0 и про-

водятся вычисления по (5.13), контролируя на каждом шаге выполнение условия монотонного убывания функции:

f (x (k +1) ) < f (x (k ) ).

Если оно выполняется, то γ0 остается неизменным; при его нарушении это число уменьшается до тех пор, пока монотонность не восстановится.

Для окончания вычислений можно использовать различные критерии, в частности,

		f (x(k) )				n		∂f	2
						= ∑				< ε.
								∂x
						i=1
								i

145

На рис. 5.8. представлена блок-схема этого алгоритма.

Начало

1

Ввод

X1H,X2H,Е,A

2

X1:=X1H

X2:=X2H

3

Y:=F(X1,X2)

4

P1:=G1(X1,X2)

P2:=G2(X1,X2)

5

R:=(P12+P22)1/2

Да

Нет

6

R<E

12

Вывод

7

XX1:=X1-A*P1

X1,X2,Y

XX2:=X2-A*P2

8

Конец

YY:=F(XX1,XX2)

Нет

Да

9

11

YY<Y

A:=A/2

Рис. 5.8

10

X1:=XX1

X2:=XX2

Y:=YY

Для простоты рассмотрен случай ЦФ двух переменных. В программе предусмотрен ввод координат начальной точки X1H, X2H, начального значения параметра γ0 = A и малого числа ε = Е.

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

146

При использовании градиентных методов (5.13) для таких функций
обычно происходит быстрый спуск со склона на дно оврага, а затем мед-
ленное	зигзагообразное		дви-	хX2
жение к точке минимума. В					2
таких	случаях гораздо		эф-
фективнее с точки зрения ус-
корения	сходимости поиска
оказывается		модификация
градиентного метода, назы-
ваемая методом сопряженных
градиентов, в частности, его
вариант,		предложенный
Флетчером и Ривсом [37]. Он					х1
состоит в выполнении сле-					X1
дующих операций:
1. Для выбранной на-
чальной точки		x(0) вычисля-			Рис. 5.9
ется антиградиент:

s (0) = −f (x (0) ),

2. На k-м шаге (k = 1, 2,…) рассматривается функция

f (x (k ) +α s (k ) )

и, используя какую-либо процедуру одномерного поиска, находится число α = α*k , при котором эта функция принимает минимальное значение; фик-

сируется получающаяся при этом точка минимума

x(k +1) = x (k) + α*k s (k), .

3. Вычисляются

f(x (k +1) ) и f (x (k +1) ) .

4.Проверяется условие прекращения поиска, например

f (x(k +1) )

∂f

2

= ∑

< ε

∂x

i

или

x(k +1) − x k )

= ∑ x(k+1)

− x(k)

2 < ε,

i

i

147

при его выполнении решение прекращается, принимая x (k+1) ≈ x*за искомое решение, иначе осуществляется переход к следующей операции.

5. Находится величина
c	/ d	,	если	k +1 I,
βk ={k	0,k		если	k +1 I,

где числа сk и dk определяются как скалярные произведения

dk = [ f (x (k) )]T [ f (x (k) )];

ck = [ f (x (k+1) )]T [ f (x (k+1) − f (x (k) )];

I ={0, n, 2n,...} – множество индексов.

6. Вычисляется s(k +1) = − f (x (k +1) ) +βk s(k) , а затем все процедура повторяется.

Реализующая эти операции укрупненная блок-схема алгоритма представлена на рис. 5.10 [37].

В ней предусмотрены: ввод начальной точки XH и допустимой погрешности E = ε, вычисление ЦФ и ее частотных производных в подпрограммах F и G, определение величин:

sp = ck ; pN = dk ; BET =βk .

Пример. Найти наименьшее значение функции Розенброка [6]: y =100 (x2 −x12 )2 +(1−x1 )2.

Заметим, что искомое решение здесь, очевидно, ymin = 0 при x1* =1, x*2 =1. Однако отыскание его численным методом в данном случае затруднено из-за того, что эта функция плохо масштабирована (1-е слагаемое в отличие от 2-го входит с большим коэффициентом), поэтому она относится к функциям овражного типа. Ее линии равного уровня показаны на рис. 5.11. Указанный особый характер этой функции при простоте аналитического выражения обусловил широкое использование ее в качестве тестовой функции при исследовании эффективности различных алгоритмов оптимизации [2, 37].

148

Начало

1

Ввод:_ N,Е,XH

2

_

_

X=XH

3

Y=F(Х)

P=G(Х)

4

K=0

5

PN=∑Pni

6

_

_

S=-P

7

*

Определение R ,миними-

зирующего F(X+R*S)

8__________

_________

__________

_

_

Х1=X+R*S

9__________

__________

_

P1=G(X1)

10

SN=∑(x1i-Xi)2

Нет

11

Да

12

15

_ _

PN=

∑

P1

2

SN<E

X*=X1

i

16

13

_

SP=∑P1i(P1i-Pi)

Y*=F(X*)

17

14

K=K+1

Вывод

X*,Y*

18

IND=K-(K/N)*N

19

конец

BET=SP/PN

Да

Нет

20

IND=0

22

PN=PN1

21

BET=0

Рис. 5.10

25 _ _

P=P1

24 _ _

X=X1

23

_

S=P1+BET S

149

Рис. 5.11

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

(x1Н = −1,5, х2Н =1,5). Поиск завершился при х1 = 0,991; х2 = 0,992.

§ 5.4. Численные методы многомерной оптимизации, без использования производных

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

150