Metody optymalizacji - teoria i wybrane

algorytmy

Michał Lewandowski

21 grudnia 2010

Spis treści

I

Algorytmy optymalizacji funkcji jednej zmiennej

3

1

Metody ustalania przedziału, w którym znajduje się minimum

4

1.1

Metoda

wyczerpującego

poszukiwania

(ang.

exhaustive se-

arch method

)

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

4

1.2

Metoda przyspieszonego poszukiwania (ang. bounding phase

method

)

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

4

2

Metody znajdowania minimum z zadaną dokładnością:

5

2.1

Metody eliminowania obszarów

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

5

Metoda dzielenia przedziału na połowę

.

.

.

.

.

.

.

.

.

.

.

.

.

.

6

Metoda złotego podziału

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

6

Metoda liczb Fibonacciego

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

7

2.2

Metody estymacji punktowej

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

8

Metoda interpolacji kwadratowej Powella

.

.

.

.

.

.

.

.

.

.

.

.

.

8

2.3

Metody oparte na gradientach

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

9

Metoda Newtona-Raphsona .

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

10

Metoda siecznych (ang. secant method) .

.

.

.

.

.

.

.

.

.

.

.

.

.

11

2.4

Porównanie metod znajdowania minimum z zadaną dokład-

nością

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

11

II

Algorytmy optymalizacji funkcji wielu zmiennych

11

METODY OPTYMALIZACJI

Michał Lewandowski

3

Metody bezpośrednich poszukiwań

12

3.1

Metoda hipersześcienna

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

13

3.2

Metoda sympleksu Neldera-Meada

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

14

3.3

Metoda kierunków sprzężonych Powella.

.

.

.

.

.

.

.

.

.

.

.

.

.

15

4

Metody gradientowe

16

4.1

Metoda Cauchy’ego najszybszego spadku

.

.

.

.

.

.

.

.

.

.

.

.

.

17

4.2

Metoda Newtona

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

18

4.3

Metoda Marquardta

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

19

III

Optymalizacja z ograniczeniami.

19

5

Teoria

19

6

Algorytmy

25

6.1

Metoda funkcji kar i barier .

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

26

2

METODY OPTYMALIZACJI

Michał Lewandowski

Moja strona internetowa: http://michallewandowski.eu

Mój adres e-mail: [e-mail]

Część I

Algorytmy optymalizacji funkcji

jednej zmiennej

Algorytmy w tym rozdziale mogą być stosowane do rozwiązywania proble-

mów minimalizacji następującego typu:

min f (x )

gdzie f (x ) jest funkcją celu a x jest zmienną rzeczywistą. Wyróżnia się dwa

podstawowe typy algorytmów:

•

Metody bezpośrednich poszukiwań (ang. direct search methods) oraz

•

Metody oparte na gradientach (ang. gradient-based methods)

Metody bezpośrednich poszukiwań wykorzystują wyłącznie wartości func-

kji celu, natomiast metody oparte na gradientach wykorzystują rownież po-

chodne pierwszego i/lub drugiego rzędu funkcji celu. Ponieważ gradienty

liczymy numerycznie, funkcja celu nie musi by c rożniczkowalna ani nawet

ciągła, aby wykorzystywać metody gradientowe. Metody opisane tutaj mogą

być stosowane również do maksymalizacji. Wystarczy po prostu sformuło-

wać i rozwiązać równoważny problem dualny min(−f (x )).

Poniższe metody zakładają, że funkcja celu jest unimodalna (ang. unimodal

function

), czyli taka, że ma tylko jedno minimum lokalne. W praktyce, dzie-

li się funkcję na przedziały, w których jest ona unimodalna i dla każdego

takiego przedziału z osobna, znajduje się minimum. My będziemy zakładać,

że szukamy minimum funkcji f na przedziale [a, b].

Minimum funkcji znajduje się w dwóch fazach:

•

Metody

ustalania

przedziału,

w

którym

znajduje

się

minimum

(ang.

bracketing methods

)

•

Metody znajdowania minimum z zadaną dokładnością:

–

Metody eliminowania obszarów (ang. region elimination methods)

3

METODY OPTYMALIZACJI

Michał Lewandowski

–

Metoda estymacji punktowej (ang. point estimation method)

–

Metody oparte na gradientach (ang. gradient based methods)

1

Metody ustalania przedziału, w którym znajduje

się minimum

W tym podrozdziale przedstawione są dwie metody:

1.1

Metoda wyczerpującego poszukiwania (ang. exhaustive se-

arch method)

Metoda ta polega na porównywaniu wartości funkcji celu dla punktów jed-

nakowo od siebie odległych. Zazwyczaj poszukiwania zaczyna się od dolne-

go ograniczenia zmiennej i w pojednynczej iteracji porównuje się wartości

trzech kolejnych punktów wykorzystując założenie unimodalności.

Algorytm

1)

Ustal

x

1

=

a

,

∆x

=

(b − a)/n

(n

jest

liczbą

punktów

pośrednich),

x

2

= x

1

+ ∆x , i x

3

= x

2

+ ∆x .

2)

Jeśli f (x

1

) ≥ f (x

2

) ≤ f (x

3

), minimum znajduje się w (x

1

, x

3

), Zakończ;

W przeciwnym przypadku x

1

= x

2

, x

2

= x

3

, x

3

= x

2

+ ∆x i przejdź do

kroku 3).

3)

Czy x

3

≤ b

? Jeśli tak, to idź do kroku 2);

Jeśli nie, to nie istnieje minimum w przedziale (a, b) lub punkt brze-

gowy (a lub b) jest puntem minimalnym.

Ostateczny

przedziałuzyskany

przy

użyciu

tego

algorytmu

to

2(b − a)/n.

Średnio

potrzeba

obliczyć

(n/2 + 2)

wartości

funkcji,

aby

uzyskać

żądaną

dokładność.

1.2

Metoda przyspieszonego poszukiwania (ang. bounding pha-

se method)

Metoda ta polega na obraniu punktu początkowego i wybraniu kierunku

poszukiwań na podstawie porównania wartości funkcji w punkcie począt-

kowym oraz dwóch wartości funkcji w punktach będących w bezpośred-

nim sąsiedztwie punktu początkowego. Później znajduje się drugi kraniec

4

METODY OPTYMALIZACJI

Michał Lewandowski

przedziału stosując wykładniczą strategię poszukiwań. Poniżej użyty jest wy-

kładnik równy 2, ale można używa˙

c jakąkolwiek inną liczbę dodatnią. Wy-

kładnik wyższy niż 1, powoduje ’przyspieszanie’ wykładnicze poszukiwań,

co zmniejsza liczbę iteracji, ale dzieje się to kosztem uzyskanej dokładno-

ści. Dla porównania w metodzie ustalania przedziału uzyskana dokładność

jest lepsza, ale ilość potrzebnych iteracji jest większa.

Algorytm

1)

Wybierz punkt początkowy x

(0)

oraz wartość ∆. Ustal k = 0.

2)

Jeśli f (x

(0)

− |

∆|) ≥ f (x

(0)

) ≥ f (x

(0)

+ |∆|), wtedy ∆ > 0;

Jeśli f (x

(0)

− |

∆|) ≤ f (x

(0)

) ≤ f (x

(0)

+ |∆|), wtedy ∆ < 0;

W pozostałych przypadkach wróć do kroku 1).

3)

Ustal x

(k+1)

= x

(k)

+ 2

k

∆.

4)

Jeśli f (x

(k+1)

) < f (x

(k)

), ustal k = k + 1 i idź do kroku 3);

W przeciwnym razie minimum znajduje się w przedziale (x

(k−1)

, x

(k+1)

),

Zakończ.

2

Metody znajdowania minimum z zadaną dokład-

nością:

2.1

Metody eliminowania obszarów

Ogólna zasada metod eliminowania obszarów jest następująca:

Rozważmy dwa punkty x

1

i x

2

, które leżą w przedziale (a, b) oraz x

1

< x

2

.

Dla problemu minimalizacji funkcji unimodalnej, można wyciągnąć nastę-

pujące wnioski:

•

Jeśli f (x

1

) > f (x

2

), to minimum nie leży w (a, x

1

)

•

Jeśli f (x

1

) < f (x

2

), to minimum nie leży w (x

2

, b

)

•

Jeśli f (x

1

) = f (x

2

), to minimum nie leży ani w (a, x

1

) ani w (x

2

, b

)

Poniżej przedstawione zostaną trzy metody:

•

Metoda dzielenia przedziału na połowę (ang. interval halving method)

•

Metoda złotego podziału (ang. golden section search)

•

Metoda liczb Fibonacciego (ang. Fibonacci search)

5

METODY OPTYMALIZACJI

Michał Lewandowski

Metoda dzielenia przedziału na połowę

Metoda

ta

polega

na

wybraniu

trzech

punktów

jednakowo

odległych

od

siebie oraz od krańców przedziału oraz wyliczeniu wartości funkcji w tych

punktach, w wyniku czego można wyeliminować połowę przedziału.

Algorytm

1)

Wybierz dolne i górne ograniczenie przedziału a i b oraz małą liczbę

ε

. Niech x

m

= (a + b)/2, L

0

= L = b − a. Wylicz f (x

m

).

2)

Ustal x

1

= a + L/4, x

2

= b − L/4. Wylicz f (x

1

) oraz f (x

2

).

3)

Jeśli f (x

1

) < f (x

m

), ustal b = x

m

; x

m

= x

1

; Przejdź do kroku 5);

W przeciwnym wypadku przejdź do kroku 4).

4)

Jeśli f (x

2

) < f (x

m

), ustal a = x

m

; x

m

= x

2

; Przejdź do kroku 5);

W przeciwnym przypadku ustal a = x

1

, b = x

2

; przejdź do kroku 5).

5)

Wylicz L = b − a. Jeśli |L| < ε , Zakończ;

W przeciwnym przypadku przejdź do kroku 2).

W każdej nowej iteracji algorytmu, potrzebne jest wyliczenie dwóch war-

tości funkcji, a przedział zmniejsza się o połowę. Po n-krotnym wyliczeniu

wartości funkcji, przedział zmniejsza się do około 0.5

n/

2

L

0

. Czyli ilość ra-

zy n, ile trzeba policzyć wartości funkcji, aby osiągnąć daną dokładność ε

można policzyć z następującego wzoru:

(0.5)

n/

2

(b − a) = ε .

Metoda złotego podziału

W metodzie tej w każdej nowej iteracji potrzeba wyliczyć tylko jedną no-

wą wartość funkcji. Idea polega na tym, że spośród dwóch punktów, które

potrzebne są, aby stosować regułę eliminowania obszarów, jeden punkt jest

zawsze poprzednim a tylko drugi punkt jest nowy. Ponadto przedział zawęża

się za każdą iteracją proporcjonalnie o tyle samo, czyli o wartość ρ, która

spełnia następującą zależność:

1 − ρ

1

=

ρ

1 − ρ

Ñ ρ

=

3 −

√

5

2

≈

0.382

Algorytm

6

METODY OPTYMALIZACJI

Michał Lewandowski

1)

Wybierz dolne i górne ograniczenie przedziału a i b oraz małą liczbę

ε

. Ustal k = 1.

2)

Ustal w

1

=

a

+ (1 − ρ)(b − a) oraz w

2

=

a

+ ρ(b − a). Wylicz f (w

1

)

oraz f (w

2

), w zależności od tego, które z nich nie było wyliczone wcze-

śniej. Wyeliminuj odpowiedni region zgodnie z regułą eliminowania

obszarów. Ustal nowe a i b.

3)

Czy |a − b| < ε ? Jeśli nie, ustal k = k + 1 i przejdź do kroku 2);

Jeśli tak, Zakończ.

W tym algorytmie, po n-krotnym wyliczeniu wartości funkcji przedział zmiej-

sza sie do (0.618)

n−

1

(a − b), zatem ilość potrzebnych wyliczeń wartości funk-

cji przy danej dokładności ε można wyliczyć z następującego wzoru:

(0.618)

n−

1

(a − b) = ε

Metoda liczb Fibonacciego

W metodzie złotego podziału proporcja zmniejszania się przedziału z iteracji

na iterację pozostaje niezmienna i wynosi 0.618. W metodzie liczb Fibonac-

ciego,

idea

jest

taka

sama

jak

w

metodzie

złotego

podziału,

z

wyjątkiem

faktu, że w metodzie liczb Fibonacciego proporcja zmniejszania się prze-

działu z iteracji na iterację zmienia się tak, aby przedział zmniejszał się w

sposób optymalny (tzn. jak najbardziej). Jeśli ρ

k

oznacza proporcje, o jaką

zmniejsza się przedział w k-tej iteracji, to w metodzie Fibonacciego zachodzi

następujący związek:

1 − ρ

k

+1

1

=

ρ

k

1 − ρ

k

(1)

Okazuje się, że wartościami ρ

k

∈

(0, 1/2], gdzie k = 1, ..., N , które minimali-

zują wyrażenie (1 − ρ

1

)(1 − ρ

2

) · · · (1 − ρ

N

) i które spełniają powyższy związek

(1), są następujące liczby:

ρ

1

=

1 −

F

N

F

N

+1

ρ

2

=

1 −

F

N−

1

F

N

...

ρ

k

=

1 −

F

N−k

+1

F

N−k

+2

...

ρ

N

=

1 −

F

1

F

2

7

METODY OPTYMALIZACJI

Michał Lewandowski

gdzie F

k

oznaczają liczby Fibonacciego. Liczby Fibonacciego mają następu-

jącą charakterystykę: F

1

= 1, F

2

= 1 i

F

k

= F

k−

1

+ F

k−

2

gdzie k = 3, 4, ....

Algorytm

1)

Wybierz dolne i górne ograniczenie przedziału a i b i ustal L = b − a.

Ustal liczbę wyliczeń wartości funkcji n. Ustal k = 1.

2)

Wylicz L

∗

k

= (1 − ρ

N−k

+1

) · · · (1 − ρ

1

)L =

F

N−k

+1

F

N

+1

. Ustal x

1

= a + L

∗

k

oraz

x

2

= b − L

∗

k

.

3)

Spośród f (x

1

) i f (x

2

) wylicz tą wartość, która nie była wcześniej po-

liczona. Wyeliminuj region według ogólnej metody eliminowania ob-

szarów. Ustal nowe a i b.

4)

Czy k = n? Jeśli nie, ustal k = k + 1 i idź do kroku 2);

W przeciwnym przypadku, Zakończ.

W tym algorytmie przedział redukuje się do

2

F

N

+1

L

do czego potrzebnych

jest n wyliczeń wartości funkcji. Zatem ilość potrzebnych wyliczeń wartości

funkcji przy danej dokładności ε można wyliczyć z następującego wzoru:

2

F

N

+1

(b − a) = ε

Jednym z minusów metody liczb Fibonacciego jest fakt, iż trzeba wylicza˙

c

kolejne liczby Fibonacciego w każdej iteracji.

2.2

Metody estymacji punktowej

W poprzednim podrozdziale omawiane były metody, które porównują war-

tości funkcji. W metodach estymacji punktowej wykorzystuje się również

wielkość różnicy wartości funkcji. Poniżej omówimy tylko metodę interpo-

lacji kwadratowej Powella.

Metoda interpolacji kwadratowej Powella

Szukamy trzy punkty x

1

< x

2

< x

3

takie, że wartości funkcji w tych punk-

tach spełniają f (x

1

) > f (x

2

) < f (x

3

). Szukamy równania wielomianu kwa-

dratowego

przechodzącego

przez

punkty

(x

1

, f

(x

1

)),

(x

2

, f

(x

2

))

i

(x

3

, f

(x

3

)).

8

METODY OPTYMALIZACJI

Michał Lewandowski

W tym celu zapisujemy ogólne równanie wielomianu kwadratowego prze-

chodzącego przez punkty x

1

i x

2

:

q

(x ) = a

0

+ a

1

(x − x

1

) + a

2

(x − x

1

)(x − x

2

)

Następnie szukamy współczynników tego wielomianu:

q

(x

1

)

=

f

(x

1

) = a

0

q

(x

2

)

=

f

(x

2

) = a

0

+ a

1

(x

2

− x

1

)

q

(x

3

)

=

f

(x

3

) = a

0

+ a

1

(x

3

− x

1

) + a

2

(x

3

− x

1

)(x

3

− x

2

)

Otrzymujemy układ trzech równań, który rozwiązujemy, aby znaleźć współ-

czynniki szukanego wielomianu:

a

0

=

f

(x

1

)

a

1

=

f

(x

2

) − f (x

1

)

x

2

− x

1

a

2

=

1

x

3

− x

2

 f

(x

3

) − f (x

1

)

x

3

− x

1

−

f

(x

2

) − f (x

1

)

x

2

− x

1



Teraz

szukamy

argumentu,

dla

którego

ten

wielomian

kwadratowy

osią-

ga minimum. Ponieważ zgodnie z naszymi założeniami a

2

>

0, minimum

znajduje się tam, gdzie pochodna równa jest zero.

q

0

(x ) = 0 Ñ a

1

+ a

2

(x − x

2

+ x − x

1

) = 0 Ñ x

∗

=

x

1

+ x

2

2

−

a

1

2a

2

W punkcie x

∗

wielomian kwadratowy q (x ) osiąga minimum. Ponieważ q (x )

jest przybliżeniem funkcji f (x ), której minimum szukamy, x

∗

jest przybli-

żeniem

wartości,

w

której

funkcja

f

(x )

osiąga

minimum.

Spośród

punk-

tów (x

1

, x

2

, x

3

, x

∗

), zatrzymujemy trzy najlepsze (innymi słowy wyrzucamy

punkt, w którym wartość funkcji f (x ) jest największa) i ponownie dokonuje-

my interpolacji kwadratowej dla tych trzech punktów i szukamy minimum

otrzymanego wielomianu. Procedura ta powtarzana jest do momentu, kiedy

osiągnięta zostanie żądana dokładność.

Alogorytm

Powella,

który

został

zarysowany

powyżej

znajduje

minimum

szybciej niż metoda złotego podziału, jeśli funkcja f (x ) nie jest skrzywiona.

Dla

mocno

niesymetrycznych

funkcji,

metoda

złotego

podziału

pozostaje

jednak lepsza.

2.3

Metody oparte na gradientach

Metody opisane dotychczas wykorzystywały tylko wartości funkcji. Meto-

dy oparte na gradientach wykorzystujá natomiast dodatkowo informację o

9

METODY OPTYMALIZACJI

Michał Lewandowski

pochodnych funkcji. Gradienty zazwyczaj oblicza się numerycznie. Używa-

jąc

metody

różnic

centralnych

(ang.

central difference method

),

liczymy

pierwszą i drugą pochodną w punkcie x

k

następująco:

f

0

(x

(k)

)

=

f

(x

(k)

+ ∆x

(k)

) − f (x

(k)

−

∆x

(k)

)

2∆x

(k)

(2)

f

00

(x

(k)

)

=

f

(x

(k)

+ ∆x

(k)

) − 2f (x

(k)

) + f (x

(k)

−

∆x

(k)

)

(∆x

(k)

)

2

(3)

Parametr ∆x

(k)

powinien być mały, na przykład może stanowić 1% wartości

x

(k)

.

Metoda Newtona-Raphsona

Załóżmy,

że

możemy

policzyć

f

(x

(k)

), f

0

(x

(k)

)

i

f

00

(x

(k)

)

w

każdym

punkcie

pomiaru funkcji x

(

k

). Możemy zdefiniować wielomian kwadratowy, którego

pierwsza i druga pochodna oraz wartść w punkcie x

(k)

są identyczne z tymi

dla funkcji f (x ). Ten wielomian ma następującą postać:

q

(x ) = f (x

(k)

) + f

0

(x

(k)

)(x − x

(k)

) +

1

2

f

00

(x

(k)

)(x − x

(k)

)

2

Zamiast minimalizować funkcję f (x ) minimalizujemy jej przybliżenie q (x ).

Warunek pierwszego rzędu na istnienie minimum jest następujący:

0 = q

0

(x ) = f

0

(x

(k)

) + f

00

(x

(k)

)(x − x

(k)

)

Nowy punkt x = x

(k+1)

spełnia zatem:

x

(k+1)

= x

(k)

−

f

0

(x

(k)

)

f

00

(x

(k)

)

Metoda Newtona-Raphsona polega na kontynuowaniu powyższej procedu-

ry do momentu, w którym pochodna f

0

(x

(k+1)

) będzie wystarczająco blisko

zera. Jeśli podstawimy g (x ) = f

0

(x ), wtedy otrzymujemy formułę do itera-

cyjnego poszukiwania rozwiązania równania g (x ) = 0:

x

(k+1)

= x

(k)

−

g

(x

(k)

)

g

0

(x

(k)

)

Metoda Newtona-Raphsona działa dobrze jeśli f

00

(x ) > 0 wszędzie. Jeśli na-

tomiast f

00

(x ) < 0 dla pewnego x , algorytm może nie zbiegać do minimum.

Jeśli zamiast analitycznych pochodnych funkcji wykorzystujemy przybliżone

pochodne (por. (2) oraz (3)), wówczas algorytm powyższy nazywamy quasi-

newtonowskim.

10

METODY OPTYMALIZACJI

Michał Lewandowski

Metoda siecznych (ang. secant method)

Jest to metoda podobna do Newtona-Raphsona. Zamiast f

00

(x

(k)

), używa się

następującego przybliżenia:

f

0

(x

(k)

) − f

0

(x

(k−1)

)

x

(k)

− x

(k−1)

Otrzymuje się wtedy przedstawiony poniżej algorytm, który nazywa się me-

todą siecznych:

x

(k+1)

= x

(k)

−

x

(k)

− x

(k−1)

f

0

(x

(k)

) − f

0

(x

(k−1)

)

f

0

(x

(k)

)

Metoda siecznych wymaga dwóch punktów startowych x

(−1)

i x

(0)

.

Tak samo, jak w przypadku metody Newtona Raphsona, metodę siecznych

można wykorzystać do znajdowania pierwiastków równania g (x ) = 0. Otrzy-

mujemy wtedy algorytm:

x

(k+1)

= x

(k)

−

x

(k)

− x

(k−1)

g

(x

(k)

) − g (x

(k−1)

)

g

(x

(k)

)

2.4

Porównanie metod znajdowania minimum z zadaną dokład-

nością

Metoda liczb Fibonacciego jest najbardziej efektywną metodą eliminacji ob-

szarów jeśli początkowy przedział, w którym leży minimum jest znany. Je-

śli nie znamy początkowego przedziału oraz pochodnych funkcji, wówczas

najlepsza powinna być metoda iterpolacji kwadratowej Powella lub metoda

quasi-newtonowska. Gdy pierwsze pochodne są dostępne, metoda siecznych

lub metoda interpolacji sześciennej

1

powinna być najbardziej efektywna. W

końcu metoda Newtona-Raphsona jest najbardziej efektywna, gdy dostępne

są informacje i pierwszych i drugich pochodnych funkcji.

1

Ta metoda nie została omówiona w rozdziale. Jest to metoda podobna do interpolacji

kwadratowej, używa jednak pierwszych pochodnych funkcji w celu zmniejszenia ilości po-

trzebnych wartości funkcji w pojedynczej iteracji.

11

METODY OPTYMALIZACJI

Michał Lewandowski

Część II

Algorytmy optymalizacji funkcji

wielu zmiennych

Dana jest funkcja wielu zmiennych: f

: R

N

Ï R

. Mówimy, że punkt

¯

x

jest

punktem

stacjonarnym,

jeśli

gradient

w

tym

punkcie

jest

zerowym

wek-

torem: ∇f ( ¯

x

) = 0. Punkt ten jest lokalnym minimum, jeśli Hesjan w tym

punkcie

∇

2

f

( ¯

x

)

jest

dodatnio

określony.

Macierz

jest

dodatnio

określona

jeśli wszystkie jej wartości własne są dodatnie: λ

i

>

0,

i

= 1, 2, ..., N

2

.

W

niniejszej

części

skryptu

omówione

będą

następujące

algorytmy

mini-

malizacji funkcji wielu zmiennych:

•

Metody bezpośrednich poszukiwań:

–

Metoda hipersześcienna (ang. evolutionary optimization)

–

Metoda sympleksowa Neldera-Meada

–

Metoda

kierunków

sprzężonych

(ang.

conjugate direction

)

Po-

wella

•

Metody gradientowe (ang. descent methods)

–

Metoda najszybszego spadku (ang. steepest descent method)

–

Metoda Newtona

–

Metoda Marquardta

–

Metoda sprzężonego gradientu Fletcher-Reevesa

–

Metody quasinewtonowskie:

∗

Metoda Davidon-Fletcher-Powella (DFP)

∗

Metoda Broyden-Fletcher-Goldfarb-Shannona (BFGS)

3

Metody bezpośrednich poszukiwań

Tak jak w przypadku optymalizacji funkcji jednej zmiennej metody bezpo-

średnich poszukiwań korzystają wyłącznie z wartości funkcji w punktach,

w przeciwieństwie do metod gradientowych, które dodatkowo wykorzystują

pochodne funkcji.

2

Macierz symetryczna, a taką jest Hesjan ma wszystkie wartości własne rzeczywiste, nie

trzeba się więc martwić o zespolone wartości własne.

12

METODY OPTYMALIZACJI

Michał Lewandowski

3.1

Metoda hipersześcienna

Algorytm potrzebuje w pojedynczej iteracji 2

N

+1 punktów, z czego 2

N

punk-

tów to są wierzchołki hipersześcianu scentrowanego na pozostałym punk-

cie. Porównuje się wartości funkcji we wszystkich tych punktach i wskazuje

się najlepszy (z najmniejszą wartością). W następnej iteracji tworzy się sze-

ścian wokół tego najlepszego punktu. Jeśli najlepszym punktem okaże się

punkt, który był środkiem danego hipersześcianu, wówczas zmniejsza się

rozmiar

sześcianu.

Proces

ten

kontynuuje

się

aż

hipersześcian

stanie

się

dostatecznie mały.

Algorytm

1)

Wybierz punkt początkowy x

(0)

oraz parametry redukcji ∆

i

dla każdej

ze zmiennych, i = 1, 2, ..., N . Wybierz parametr zakończenia ε . Ustal

¯

x

= x

(0)

2)

Jeśli ||∆|| < ε , Zakończ;

W przeciwnym przypadku stwórz 2

N

punktów poprzez dodanie i od-

jęcie ∆

i

/

2 do/od każdej zmiennej w punkcie ¯

x

.

3)

Oblicz wartość funkcji dla wszystkich (2

N

+ 1) punktów. Znajdź punkt

z najmniejszą wartością. Ustal ten punkt jako ¯

x

.

4)

Jeśli

¯

x

=

x

(0)

,

zredukuj

parametry

redukcji

∆

i

=

∆

i

/

2

i

przejdź

do

kroku 2);

W przeciwnym przypadku ustal x

(

0) = ¯

x

i przejdź do kroku 2)

W powyższym algorytmie w każdej iteracji trzeba policzyć maksymalnie 2

N

wartości fukcji. Czyli ilość potrzebnych ewaluacji funkcji wzrasta wykładni-

czo wraz z N .

Jeśli

algorytm

znajdzie

dokładne

minimum

funkcji

w

którejś

iteracji,

nie

zatrzymuje się automatycznie. Wartość ||∆|| musi spaść poniżej ε , aby algo-

rytm zakończył działanie.

Ustalenie dużego parametru redukcji ∆

i

jest dobre, ale wtedy algorytm mo-

że potrzebować wielu iteracji a zatem i wielu ewaluacji funkcji. Z drugiej

strony ustalenie małego parametru redukcji może prowadzić do zbyt wcze-

snej zbieżności algorytmu do suboptymalnego punktu, szczególnie w przy-

padku funkcji bardzo nieliniowych. Redukcja parametru redukcji nie musi

być dwukrotna tak jak w algorytmie podanym powyżej (patrz krok 4)). Dla

13

METODY OPTYMALIZACJI

Michał Lewandowski

lepszej zbieżności algorytmu rekomenduje się redukcję mniejszą niż dwu-

krotną

3

.

3.2

Metoda sympleksu Neldera-Meada

Sympleks wymaga dużo mniejszej ilości punktów niż hypersześcian, co sta-

je się szcególnie widoczne dla wielu wymiarów. Sympleks w N -wymiarowej

przestrzeni

ma

N

+ 1

wierzchołków.

Jest

to

minimalna

ilość

wierzchoł-

ków umożliwiająca poszukiwanie we wszystkich możliwych kierunkach N -

wymiarowej przestrzeni. Ważne jest jednak, aby sympleks

4

nie rozpinał fi-

gury o zerowej objętości w N -wymiarowej przestrzeni - czyli na przykład

dla funkcji dwuwymiarowej trzy punkty sympleksu nie mogą leżeć na jed-

nej linii, a w przypadku funkcji trójwymiarowej cztery punkty sympleksu

nie mogą leżeć na jednej płaszczyźnie.

Działanie metody najlepiej widać w pseudokodzie algorytmu:

Algorytm

1)

Wybierz γ > 1, β ∈ (0, 1) oraz parametr zakończenia ε . Stwórz sym-

pleks początkowy.

2)

Znajdź najgorszy punkt x

w

, najlepszy punkt x

b

oraz drugi w kolejności

najgorszych x

m

. Oblicz punkt, wzgędem którego będziemy odbijać x

w

:

x

c

=

1

N

N

+1

X

i

=1,i6=h

x

i

.

3)

Oblicz punkt odbicia x

r

= 2x

c

− x

w

. Ustal x

new

= x

r

.

Jeśli f (x

r

) < f (x

b

), ustal x

new

= (1 + γ )x

c

− γx

w

(ekspansja);

Jeśli f (x

r

) ≥ f (x

w

), ustal x

new

= (1 − β)x

c

+ βx

w

(kontrakcja);

Jeśli f (x

m

) < f (x

r

) < f (x

w

), ustal x

new

= (1 + β)x

c

− βx

w

(kontrakcja).

Policz f (x

new

) i wymień x

w

na x

new

.

4)

Jeśli

n

P

N

+1

i

=1

(f (x

i

)−f (x

c

))

2

N

+1

o

1/2

≤ ε

, Zakończ;

W przeciwnym przypadku idź do kroku 2).

Jednym ze sposobów stworzenia sympleksu początkowego do pierwszego

kroku algorytmu jest wybranie punktu bazowego x

(0)

oraz liczby C . Wów-

3

W kroku 4): ∆

i

= ∆

i

/p

, gdzie p ∈ (1, 2).

4

Chodzi tutaj o sympleks początkowy. O następne nie trzeba się martwić, co będzie wi-

doczne w algorytmie.

14

METODY OPTYMALIZACJI

Michał Lewandowski

czas N + 1 punktów sympleksu to x

(0)

oraz dla i, j = 1, 2, ..., N :

x

(i)

j

=

(

x

(0)

j

+ C

if j = i

x

(0)

j

+ C ∆

otherwise,

where ∆ =

(

0.25

if N = 3

√

N

+1−2

N−

3

otherwise.

Dla zapewnienia dobrej zbieżności algorytmu, sugeruje się ustalenie γ ≈ 2

i |β| ≈ 0.5.

3.3

Metoda kierunków sprzężonych Powella.

Metoda kierunków sprzężonych Powella jest chyba najbardziej popularną

metodą bezpośrednich poszukiwań. Wykorzystuje ona historię poprzednich

rozwiązań, aby swtorzyć nowe kierunki poszukiwań. Idea jest prosta: trzeba

utworzyć N liniowo niezależnych kierunków poszukiwań i dokonać sekwen-

cyjnie serię poszukiwań wzdłuż tych kierunków, startując za każdym razem

z poprzednio znalezionego punktu.

Algorytm

1)

Wybierz punkt początkowy x

(0)

oraz zbiór N liniowo niezależnych kie-

runków; najlepiej s

(i)

= e

(i)

, dla i = 1, 2, ..., N , gdzie e

(i)

oznacza i-ty

wektor bazowy bazy kanonicznej.

2)

Szukaj minimum startując z punktu początkowego wzdłuż kierunku

s

(1)

. Startując z nowo znalezionego punktu (oznacz go jako y

(1)

), szukaj

minimum wzdłuż kierunku s

(2)

. Kontynuuj szukanie wzdłuż kolejnych

kierunków aż do kierunku s

(N )

. Następnie ponownie szukaj minimum

wdłuż kierunku s

(1)

. Punkt, który został ostatecznie znaleziony oznacz

jako y

(2)

.

3)

Oblicz d = y

(2)

− y

(1)

. Jest to kierunek sprzężony do s

(1)

.

4)

Jeśli ||d|| jest małe lub kierunki s

(1)

, s

(2)

, ..., s

(N −1)

, d

są liniowo zależne,

Zakończ

;

W przeciwnym przypadku ustal s

(j )

= s

(j −1)

dla wszystkich j = N , N −

1, ..., 2. Ustal s

(1)

= d/||d|| i idź do kroku 2).

15

METODY OPTYMALIZACJI

Michał Lewandowski

4

Metody gradientowe

Gradient w punkcie x

(t )

możemy przybliżyć numerycznie za pomocą nastę-

pującej formuły:

∇f

(x

(t )

) =









∂f

(x

(t )

)

∂x

1

∂f

(x

(t )

)

∂x

2

...

∂f

(x

(t )

)

∂x

N









,

gdzie

∂f

(x

(t )

)

∂x

i

=

f

(x

(t )

i

+ ∆x

(t )

i

) − f (x

(t )

i

−

∆x

(t )

i

)

2∆x

(t )

i

Hesjan w punkcie x

(t )

natomiast liczymy następująco:

∇

2

f

(x

(t )

=











∂

2

f

(x

(t )

)

∂x

2

1

∂

2

f

(x

(t )

)

∂x

1

∂x

2

...

∂

2

f

(x

(t )

)

∂x

1

∂x

N

∂

2

f

(x

(t )

)

∂x

2

∂x

1

∂

2

f

(x

(t )

)

∂x

2

2

...

∂

2

f

(x

(t )

)

∂x

2

∂x

N

...

...

...

...

∂

2

f

(x

(t )

)

∂x

N

∂x

1

∂

2

f

(x

(t )

)

∂x

N

∂x

2

...

∂

2

f

(x

(t )

)

∂x

2

N











gdzie:

∂

2

f

(x

(t )

)

∂x

2

i

=

f

(x

(t )

i

+ ∆x

(t )

i

) − 2f (x

(t )

) + f (x

(t )

i

−

∆x

(t )

i

)

(∆x

(t )

i

)

2

∂

2

f

(x

(t )

)

∂x

i

∂x

j

=

∂f

(x

(t )

i

+∆x

(t )

i

)

∂x

j

−

∂f

(x

(t )

i

−

∆x

(t )

i

)

∂x

j

2∆x

(t )

i

Pochodne cząstkowe w ostatnim wyrażeniu powyżej są z kolei liczone tak,

jak składowe gradientu, tylko że w innym punkcie. Wyrażenie x

(t )

i

+ ∆x

(t )

i

reprezentuje wektor (x

(t )

1

, ..., x

(t )

i

+ ∆x

(t )

i

, ..., x

(t )

N

)

T

Żeby policzyć gradient potrzebnych jest 2N różnych wartości funkcji. Ażeby

policzyć Hesjan potrzebnych jest 3N + 4

N

2

!

= 2N

2

+ N .

Ponieważ gradient jest kierunkiem najszybszego wzrostu, minus gradient

jest kierunkiem najszybszego spadku funkcji. Kierunek poszukiwań (ang.

search direction

) d

(t )

jest kierunkiem spadku w punkcie x

(t )

, jeśli w otocze-

niu tego punktu spełniony jest następujący warunek:

∇f

(x

(t )

) · d

(t )

≤

0

Oznacza to, że cosinus kąta między gradientem i kierunkiem poszukiwań

jest

większy

niż

90

0

.

Kierunek

d

(t )

jest

kierunkiem

spadku,

ponieważ

w

wyniku rozwinięcia f wokół x

(t )

otrzymujemy:

f

(x

(t +1)

) = f (x

(t )

+ αd

(t )

) = f (x

(t )

) + α∇f (x

(t )

) · d

(t )

.

16

METODY OPTYMALIZACJI

Michał Lewandowski

Im niższa ujemna wartość ∇f (x

(t )

) · d

(t )

tym większy spadek funkcji w kie-

runku d

(t )

.

W metodach gradientowych często w ramach pojedynczej iteracji dokonuje

się poszukiwań w danym kierunku. Jest to optymalizacja jednowymiarowa.

Najpierw zapisujemy reprezentatywny punkt wzdłuż kierunku s

(t )

jako:

x

(k+1)

= x

(k)

+ α

(k)

s

(k)

,

gdzie α

(k)

jest długością kroku. Ponieważ x

(k)

i s

(k)

są znane, punkt x

(k+1)

można

zapisać

tylko

jedną

zmienną.

Można

więc

wykonać

minimalizację

jednowymiarową, aby otrzymać nowy punkt x

(k+1)

. Następnie poszukiwania

kontynuuje się wzdłuż nowego kierunku s

(k+1)

i tak dalej aż do momentu

znalezienia lokalnego minimum. Jeśli metoda gradientowa jest użyta do po-

szukiwań jednowymiarowych wzdłuż kierunku, minimum znajdujemy po-

przez

zróżniczkowanie

wyrażenia

f

(x

(t +1)

)

=

f

(x

(t )

+ αs

(k)

)

względem

α

i

przyrównaniem do zera:

∇f

(x

(k+1)

) · s

(k)

= 0.

W ten sposób znajdujemy nowy punkt x

(k+1)

. Okazuje się że kąt pomiędzy

kierunkiem poszukiwań w k-tej iteracji i kierunkiem najszybszego spadku

w nowym punkcie −∇f (x

(k+1)

) jest równy 90

0

.

4.1

Metoda Cauchy’ego najszybszego spadku

Kierunek poszukiwań w metodzie Cauchy’ego jest kierunkiem najszybsze-

go spadku:

s

(k)

= −∇f (x

(k)

).

Algorytm ten gwarantuje poprawę, to jest spadek wartości funkcji, w każdej

iteracji. Metoda najszybszego spadku działa dobrze, gdy x

(0)

jest daleko od

minimum x

∗

. Jeśli bieżący punkt jest blisko minimum, zmiana gradientu jest

mała, wobec następny punkt powstały w wyniku poszukiwania w jednym

kierunku jest blisko punktu bieżącego.

Algorytm

1)

Wybierz maksymalną liczbę iteracji M , punkt początkowy oraz dwa

parametry zakończenia ε

1

i ε

2

i ustal k = 0.

2)

Oblicz ∇f (x

(k)

)

17

METODY OPTYMALIZACJI

Michał Lewandowski

3)

Jeśli ||∇f (x

(k)

)|| ≤ ε

1

, Zakończ;

Jeśli k ≥ M ; Zakończ;

W przeciwnym razie idź do kroku 4).

4)

Wykonaj

poszukiwanie

wzdłuż

kierunku,

żeby

znaleźć

α

(k)

tak,

aby

f

(x

(k+1)

)

=

f

(x

(k)

+ α

(k)

s

(k)

) było minimalne. Jednym z kryteriów za-

kończenia jest |∇f (x

(k+1)

) · ∇f (x

(k)

)| ≤ ε

2

.

5)

Jeśli

||x

(k+1)

−x

(k)

||

||x

(k)

||

≤ ε

1

, Zakończ;

W przeciwnym przypadku ustal k = k + 1 i idź do kroku 2).

4.2

Metoda Newtona

Rozwinięcie w szereg Taylora drugiego rzędu funkcji f

wokół punktu x

(t )

przyjmuje następującą formę:

f

(x

(k+1)

) = f (x

(k)

) + ∇f (x

(k)

)

T

(x − x

(k)

) +

1

2

(x − x

(k)

)

T

∇

2

f

(x

(k)

)(x − x

(k)

)

Jeśli policzymy warunek pierwszego rzędu na maksimum lokalne tej funk-

cji, otrzymujemy:

∇f

(x

(k)

) + ∇

2

f

(x

(k)

)(x − x

(k)

) = 0

Podstawiając x

(k+1)

= x , otrzymujemy:

x

(k+1)

= x

(k)

−

h

∇

2

f

(x

(k)

)

i

−

1

∇f

(x

(k)

)

Kierunek poszukiwań w metodzie Newtona jest zatem dany wyrażeniem:

s

(k)

= −

h

∇

2

f

(x

(k)

)

i

−

1

∇f

(x

(k)

)

Jeśli macierz

h

∇

2

f

(x

(k)

)

i

−

1

jest półdodatnio określona, kierunek s

(k)

jest kie-

runkiem spadku. Warunek drugiego rzędu optymalizacji mówi, że macierz

∇

2

f

(x

∗

) jest dodatnio określona dla minimum lokalnego. Można zatem za-

łożyć,

że

macierz

∇

2

f

(x )

jest

dodatnio

określona

w

otoczeniu

minimum.

Metoda Newtona jest więc dobra, kiedy punkt początkowy jest blisko mini-

mum.

Algorytm jest bardzo podobny do metody najszybszego spadku. Poszukiwa-

nia prowadzone są jednak w innym kierunku s

(k)

= −

h

∇

2

f

(x

(k)

)

i

−

1

∇f

(x

(k)

).

Możliwy warunek zakończenia optymalizacji wzdłuż kierunku wygląda na-

stępująco:






∇f

(x

(k+1)

) ·

h

∇

2

f

(x

(k)

)

i

−

1

∇f

(x

(k)

)






≤ ε

2

.

18

METODY OPTYMALIZACJI

Michał Lewandowski

4.3

Metoda Marquardta

Metoda Cauchy’ego działa dobrze, gdy punkt początkowy jest daleko od mi-

nimum, podczas, gdy metoda Newtona działa dobrze, gdy punkt początkowy

jest blisko minimum. W metodzie Marquardta, metodę Cauchy’ego stosuje

się na początku by następnie zaadoptować metodę Newtona.

Algorytm

1)

Wybierz maksymalną liczbę iteracji M , punkt początkowy oraz para-

metr zakończenia ε . Ustal k = 0 oraz λ

(0)

= 10

4

(duża liczba).

2)

Oblicz ∇f (x

(k)

)

3)

Jeśli ||∇f (x

(k)

)|| ≤ ε , Zakończ;

Jeśli k ≥ M ; Zakończ;

W przeciwnym razie idź do kroku 4).

4)

Oblicz x

(k+1)

= x

(k)

−

h

∇

2

f

(x

(k)

) + λ

(k)

I

i

−

1

∇f

(x

(k)

)

5)

Jeśli f (x

(k+1)

) < f (x

(k)

), idź do kroku 6);

W przeciwnym przypadku idź do kroku 7).

6)

Ustal λ

(k+1)

=

1

2

λ

(k)

, k = k + 1 i idź do kroku 2).

7)

Ustal λ

(k+1)

= 2λ

(k)

i idź do kroku 4).

Algorytm może być szybszy, jeśli dodatkowo będziemy dokonywać optyma-

lizacji wzdłuż kierunku w każdej iteracji.

Część III

Optymalizacja z ograniczeniami.

5

Teoria

Nieniejszy rozdziałma służyć intuicyjnemu przedstawieniu idei twierdzenia

Lagrangeá i Kuhn-Tuckera.

Dane jest następujące zadanie optymalizacyjne:

max

x∈D⊂R

N

f

(x )

19

METODY OPTYMALIZACJI

Michał Lewandowski

Zbiór D jest zbiorem punktów dopuszczalnych. Zbiór ten będziemy przed-

stawiać za pomocą ograniczeń w postaci równości h(x ) = c oraz ograniczeń

w postaci nierówności g (x ) ≤ c. W przypadku ograniczeń w postaci nierów-

ności, mówimy, że dane ograniczenie jest aktywne bądź napięte w punkcie

dopuszczalnym x jeśli zachodzi g (x ) = c. W przeciwnym przypadku, tj. gdy

g

(x ) < c mówimy, że ograniczenie jest nieaktywne lub luźne w punkcie do-

puszczalnym x . Oczywiście dla ograniczeń w postaci równości dla punktu

dopuszczalnego ograniczenie z definicji musi być napięte. Intuicję dotyczącą

zagadnień optymalizacyjnych najlepiej wyrobić sobie graficznie dla przypad-

ku, gdy x ∈ R

2

.

Rysunki poniżej pokazują, że w przypadku ograniczeń w postaci równości w

punkcie optymalnym

5

nachylenie funkcji celu i ograniczenia powinno być

równe:

∇h(x*)

x

2

x*

x'

∇h(x')

∇f(x')

x

1

A

h(x)=c

∇f(x*)

∇f(x*)=λ∇h(x*), λ>0

Równe nachylenie i mnożnik dodatni

∇h(x*)

x

2

x*

x'

∇h(x')

∇f(x')

x

1

A

h(x)=c

∇f(x*)

∇f(x*)=λ∇h(x*), λ<0

Równe nachylenie i mnożnik ujemny

5

Zakładając różniczkowalność.

20

METODY OPTYMALIZACJI

Michał Lewandowski

Możemy to zapisać w następujący sposób:

∇f

(x

∗

) = λ∇h(x

∗

)

(4)

gdzie symbol λ oznacza mnożnik Lagrange’a.

Tak więc problem optymalizacji z ograniczeniem max

x

f

(x ),

p.w.

h

(x ) =

c

można

sprowadzić

do

problemu

optymalizacji

bez

ograniczeń

funkcji

L

(x , λ)

= f (x ) − λ(h(x ) − c), zwanej funkcją Lagrange’a, bowiem warunki

pierwszego rzędu optymalizacji tej funkcji będą dokładnie równe warun-

kom (4) oraz ograniczeniu h(x

∗

) = c:

∇f

(x

∗

) − λ∇h(x

∗

)

=

0

h

(x

∗

) − c

=

0

Z powyższych rysunków wynika również, że w przypadku ograniczeń w po-

staci równości mnożnik Lagrange’a może być zarówno dodatni, jak i ujemny.

Jedyna wymagana rzecz jest taka, że musi on być zdefiniowany. Aby mnoż-

nik można było zdefiniować musi być spełniony następujący warunek:

∇h

(x

∗

) 6= 0

(5)

Gradient ograniczenia w punkcie optymalnym musi być różny od wektora

zerowego. W naszym przykładzie oznacza to, że:

"

∂h

(x

∗

)

∂x

1

∂h

(x

∗

)

∂x

2

#

6

=

"

0

0

#

Dla problemu optymalizacji z ograniczeniami w postaci nierówności g (x ) ≤

c

problem jest podobny, lecz nieco bardziej złożony. Jeśli optimum jest w

punkcie, w którym dane ograniczenie nie jest aktywne, to warunek pierw-

szego rzędu jest taki sam jak, gdyby ograniczenia nie było tj.

∇f

(x

∗

) = 0

(6)

Warunek pierwszego rzędu funkcji Lagrange’a wynosi zaś:

∇f

(x

∗

) − λg (x

∗

) = 0

Aby te warunki się zgadzały, mnożnik Lagrange’a w tym przypadku musi

wynosić zero.

Jeśli zaś optimum jest w punkcie, w którym dane ograniczenie jest aktywne,

wówczas mamy do czynienia z sytuacją jak na rysunku poniżej:

21

METODY OPTYMALIZACJI

Michał Lewandowski

∇f(x*)

x

2

x

1

x*

x'

∇g(x')

∇f(x')

A

g(x)≤c

∇g(x*)

∇f(x*)=λ∇g(x*), λ>0

W tym przypadku mnożnik Lagrange’a musi być dodatni, ponieważ w prze-

ciwnym razie moglibyśmy przesunąć punkt x

∗

do wewnątrz zbioru dopusz-

czalnego podnosząc wartość funkcji celu f

- co przeczyłoby optymalności

punktu x

∗

:

∇g(x*)

∇f(x*)

x

2

x

1

Dlaczego mnożnik Lagrange'a

musi być dodatni

g(x)=b

g(x)>b

g(x)<b

Zatem mamy dwie sytuacje możliwe w punkcie optymalnym:

•

ograniczenie aktywne, mnożnik dodatni: g (x ) = c i λ > 0

•

ograniczenie nieaktywne, mnożnik zerowy: g (x ) < 0 i λ = 0

22

METODY OPTYMALIZACJI

Michał Lewandowski

Te dwa przypadki można zapisać jednym warunkiem zwanym "complemen-

tary slackness":

λ

[g (x

∗

) − c] = 0

Analiza będzie podobna jeśli ograniczenia w postaci ≤ zamienimy na ogra-

niczenia w postaci ≥. Wówczas trzeba jednak w prezentacji graficznej zmo-

dyfikować kierunki gardientów ograniczeń:

x

2

x

1

{x|g

2

(x)=b

2

}

∇g

1

(x)

∇g

4

(x)

∇g

3

(x)

∇g

2

(x)

g

1

(x)≥b

1

g

2

(x)≥b

2

g

3

(x)≥b

3

g

4

(x)≥b

4

Ograniczenia w postaci ≥

x

2

x

1

{x|g

2

(x)=b

2

}

∇g

1

(x)

∇g

4

(x)

∇g

3

(x)

∇g

2

(x)

g

1

(x)≤b

1

g

2

(x)≤b

2

g

3

(x)≤b

3

g

4

(x)≤b

4

Ograniczenia w postaci ≤

W przypadku więcej niż jednego ograniczenia, zarówno jeśli chodzi o ogra-

niczenia w postaci równości, które w punkcie optymalnym muszą być ak-

23

METODY OPTYMALIZACJI

Michał Lewandowski

tywne

z

definicji

6

,

jak

i

ograniczenia

w

postaci

nierówności,

ale

które

w

punkcie optymalnym są aktywne, mamy do czynienia z sytuacją, jak na ry-

sunku poniżej:

x

2

x

1

∇g

1

(x*)

∇g

2

(x*)

g

1

(x)=c

1

g

1

(x)≤c

1

g

2

(x)=c

2

g

2

(x)≤c

2

∇f(x*)

∇f(x*)=λ

1

∇g

1

(x*)+λ

2

∇g

2

(x*)

Kombinacja liniowa gradientów ograniczeń aktywnych

Tym razem gradient funkcji celu musi być kombinacją liniową gradientów

ograniczeń

aktywnych.

W

punkcie

x

∗

oba

ograniczenia

g

1

(x )

≤ c

1

oraz

g

2

(x ) ≤ c

2

są aktywne, a zatem w tym punkcie gradient funkcji celu jest

kombinacją liniową gradientów obu ograniczeń. Mnożniki Lagrange’a od-

powiadające poszczególnym ograniczeniom występują tutaj w roli wag kom-

binacji liniowej. Oczywiście, aby gradient funckji celu dało się zapisać jako

kombinacjłe liniową gradientów ograniczeń aktywnych, te gradienty ogra-

niczeń aktywnych powinny być liniowo niezależne. Inaczej możemy mieć

do czynienia z sytuacją, jak na rysunku poniżej:

6

Punkt optymalny musi być punktem dopuszczalnym.

24

METODY OPTYMALIZACJI

Michał Lewandowski

∇f(x*)

∇g

1

(x*)

∇g

2

(x*)

Dlaczego potrzebne jest constraint qualification?

x

2

x

1

Jedyny

punkt

dopuszczalny

to

punkt

x

∗

,

zatem

jest

on

zarazem

punktem

optymalnym. Jednakże nie da się zapisać ∇f (x

∗

) jako liniową kombinację

gradientów ∇g

1

(x

∗

) i ∇g

2

(x

∗

). Liniowa niezależność gradientów wszystkich

ograniczeń w postaci równości oraz aktywnych ograniczeń w postaci nie-

równości jest zwana warunkiem "constraint qualification".

Warunek konieczny nie wystarczajacy

Nie mozna pozbywac sie nieaktywnych ograniczen

Dlaczego x=0 nie mozna zastapic x<=0 i -x<=0.

Cofnij sie do h(x)=c i wytlumacz ten punkcik A na rysuneczkach

6

Algorytmy

Dane są funkcje f : R

N

Ï R

, g

j

: R

N

Ï R

, gdzie j = 1, ..., J oraz h

k

: R

N

Ï R

,

gdzie k = 1, ..., K. Ogólny problem optymalizacyjny w niniejszej części jest

następujący:

min

x

f

(x )

przy warunkach:

g

j

(x ) ≥ 0,

j

= 1, ..., J ;

h

k

(x ) = 0,

k

= 1, ..., K;

Funkcja Lagrange’a dla powyższego problemu ma postać:

L

(x , u, v ) = f (x ) −

J

X

j

=1

u

j

g

j

(x ) −

K

X

k

=1

v

k

h

k

(x )

25

METODY OPTYMALIZACJI

Michał Lewandowski

Warunki Kuhn-Tuckera zapisujemy następująco:

∇f

(x ) −

J

X

j

=1

u

j

∇g

j

(x ) −

K

X

k

=1

v

k

∇h

k

(x ) = 0

g

j

(x ) ≥ 0,

j

= 1, ..., J ;

h

k

(x ) = 0,

k

= 1, ..., K;

u

j

g

j

(x ) = 0,

j

= 1, ..., J ;

u

j

≥

0,

j

= 1, ..., J .

Możemy teraz zapisać twierdzenie Kuhn-Tuckera o warunku koniecznym

istnienia optimum.

Twierdzenie.

Dla problemu zadanego powyżej, niech f , g

j

oraz h

k

bę-

dą różniczkowalne a punkt x

∗

punktem dozwolonym. Niech I

(x

∗

) = {j

:

g

j

(x

∗

) = 0} oznacza zbiór aktywnych ograniczeń w postaci nierówności.

Niech ∇g

j

(x

∗

) dla j ∈ I oraz ∇h

k

(x

∗

), dla k = 1, ..., K będą liniowo nieza-

leżne ("constraint qualification"). Jeśli x

∗

jest optymalnym rozwiązaniem

problemu, wówczas istnieje wektor mnożników Lagrange’a

(u

∗

, v

∗

), taki że

(x

∗

, u

∗

, v

∗

) spełnia warunki Kuhn-Tuckera.

Warunek "constraint qualification" oznacza, że gradienty wszystkich ogra-

niczeń aktywnych w dozwolonym punkcie są liniowo niezależne

7

.

Zatem każdy optymalny punkt musi spełniać warunki Kuhn-Tuckera, ale

nie każdy punkt, który spełnia warunki Kuhn-Tuckera jest optymalny. Jeśli

chodzi o punkty, które nie spełniają "constraint qualification", nic się nie da

o nich powiedzieć - mogą być optymalne i mogą być nieoptymalne.

6.1

Metoda funkcji kar i barier

Algorytm funkcji kar i barier polega na sprowadzeniu problemu optyma-

lizacji z ograniczeniami do serii problemów optymalizacji bez ograniczeń.

Dla

różnych

parametrów

funkcji

kar

i

barier

wykonuje

się

wielokrotnie

optynalizację funkcji bez ograniczeń, która dana jest nastę pują cym wzo-

rem:

P

(x , R) = f (x ) + Ω(R, g (x ), h(x ))

gdzie R jest zbiorem parametrów kary a Ω jest funkcją kary, która fawory-

zuje selekcję punktów dopuszczalnych nad punkty niedopuszczalne. Używa

7

Ograniczenia w postaci równości są aktywne z definicji dla punktu dopuszczalnego, a

spośród pozostałych ograniczeń - tych w postaci nierówności - aktywne są tylko te ze zbioru

I

(x

∗

).

26

METODY OPTYMALIZACJI

Michał Lewandowski

się różnych funkcji kary dla ograniczeń w postaci równości i dla ograniczeń

w postaci nierówności. Dodatkowo niektóre funkcje Ω nakładają karę tylko

na punkty niedozwolone (wówczas nazywają się funkcjami kary lub funkcja-

mi kary zewnętrznej) a niektóre funkcje nie potrafią w ogóle radzić sobie z

punktami niedopuszczalnymi i nakładają karę na punkty dozwolone, jeśli są

blisko granicy ograniczenia (nazywają się wówczas funkcją bariery lub we-

wnętrzną funkcją kary). Poniżej przedstawiona będzie jedna wersja metody

kar i barier, w której będziemy zajmować się wyłą cznie ograniczeniami w

postaci nierówności i funkcja kary bedzie nastę pują cej postaci:

Ω = Rhg (x )i

2

gdzie hg (x )i =

(

g

(x ),

jeśli g (x ) < 0

0,

jeśli g (x ) ≥ 0

Na początku optymalizuje się f (x ) + Ω dla małej wartości R, następnie sukce-

sywnie zwiększa się wartość R i znowu optymalizuje metodami optymalizacji

bez ograniczeń dla nowych wartości R. Poniżej znajduje się algorytm:

Algorytm

1)

Wybierz parametry zakończenia ε

1

, ε

2

, początkowy punkt x

(0)

, począt-

kowy parametr funkcji kary R

(0)

oraz parametr zwiększania parame-

tru funkcji kary c > 1. Ustal t = 0.

2)

Utwórz P (x

(t )

, R

(t )

) = f (x

(t )

) + Rhg (x

(t )

)i

2

.

3)

Startując z x

(t )

, znajdź x

(t +1)

tak, że P (x

(t +1)

, R

(t )

) jest minimalne dla

ustalonej wartości R

(t )

- wykorzystaj jedną z metod minimalizacji funk-

cji wielu zmiennych bez ograniczeń. Do zatrzymania poszukiwań użyj

ε

1

.

4)

Jeśli |P (x

(t +1)

, R

(t )

) − P (x

(t )

, R

(t −1)

)| ≤ ε

2

, ustal x

T

= x

(t +1)

i Zakończ.

W przeciwnym przypadku idź do kroku 5).

5)

Wybierz R

(t +1)

= cR

(t )

. Ustal t = t + 1 i idź do kroku 2).

27