Strona główna

Zajęcia : Optymalizacja funkcji wielu zmiennej


Pobieranie 48.46 Kb.
Data19.06.2016
Rozmiar48.46 Kb.

Zajęcia 3.: Optymalizacja funkcji wielu zmiennej.


Problem optymalizacji

pw. n

f(x) – funkcja unimodalna

3.1 Metody bezpośrednich poszukiwań1

3.1.1 Metoda pełzającego sympleksu Nelder’a-Mead’a


Definicja

Simpleksem nazywamy figurę geometryczna składającą się z N+1 punktów w przestrzeni N-wymiarowej, np. sympleksem w przestrzeni R2 jest trójkąt, zaś w R3 czworościan. Simpleks o równej długości boków nazywany jest simpleksem regularnym.

Dane:


f(x) – funkcja N zmiennych

γ = (1, ∞) – czynnik ekspansji/rozciągnięcia (ang. Expansion phase)

β = (0, 1) – czynnik kontrakcji/ściśnięcia (ang. Contraction phase)

{x1, x2, …, xN, xN+1} – początkowy sympleks, czyli zbiór N+1 wierzchołków figury

ε – parametr stopu


  1. Ze zbioru punktów simpleksu znajdź:

    1. Najgorszy punkt i oznacz go jako xh

    2. Drugi z najgorszych i oznacz go jako xg

    3. Najlepszy punkt i oznacz go jako xl

Policz środek ciężkości/symetrii dla punktów simpleksu z wyłączeniem punktu xh

Xc = ∑xi / N (dla i ≠ h)

  1. Oblicz punkt odbicia2 xr punktu xh:

xr = xc + (xcxh) = 2xc - xh

Ustal xnew jako xr, tj. xnew = xr (ani ekspansja ani kontrakcja)

Jeżeli f(xr) < f(xl), to xnew = xc + γ (xcxh) = (1 + γ) xc - γ xh (ekspansja)

Jeżeli f(xr) > f(xh), to xnew = xc – β (xcxh) = (1 – β) xc + β xh (kontrakcja)

Jeżeli f(xg) < f(xr) < f(xh), to xnew = xc + β (xcxh) = (1 + β) xc - β xh (kontrakcja)

Zastąp w simpleksie, xh przez xnew.



  1. Jeżeli { ∑ (f(xi) – f(xc))2 / (N +1) }1/2 < ε, zakończ algorytm, wpp.  pkt. 1

Uwagi:

- Metoda Nelder’a – Mead’a jest stosowana najczęściej do funkcji nieróżcznikowalnych albo nieciągłych. W przypadku funkcji różniczkowalnych i dostępności informacjo o pochodnych, zaleca się stosowanie metod gradientowych.

- W przypadku gdyby sympleks nabrał nieregularnych kształtów, tj. był za bardzo rozciągnięty w jedną tylko stronę zaleca się zastąpienie takiego sympleksu bardziej regularnym. Procedura taka nazywa się periodyczną odnową.

- duże wartości γ i 1 / β sprzyjają szybszemu zbliżaniu się do optimum w pierszych krokach metody, ale utrudniają jej zbieganie pod koniec algorytmu. Z drugiej strony małe wartości γ i 1 / β mogą wymagać dużej liczby obliczeń wartości funkcji celu.

- zalecane wartości parametrów metody to: γ = 2 i β = 0,5

3.1.2 Metoda kierunków sprzężonych Powella


Własność równoległych podprzestrzeni

Mamy daną funkcję kwadratową dwóch zmiennych f(x) = a + bTx + ½ xTCx oraz dwa punkty x1x2 oraz kierunek d.

Jeżeli y1 jest rozwiązaniem problemu minimalizacji w kierunku min f(x1 + λ d) oraz y2 jest rozwiązaniem problemu minimalizacji w kierunku min f(x2 + λ d), to kierunek (y2y1) jest sprzężony do d, innymi słowy (y2y1)TC d = 0.

Uogólniona własność równoległych podprzestrzeni

Załóżmy, że y1 i y2 został znaleziony po m przeszukiwaniach w kierunku zaczynających się w punktach odpowiednio: x1 i x2. Wektor (y1y2) jest sprzężony do wszystkich m kierunków poszukiwań.


Dane:

f(x) – funkcja N zmiennych



x0punkt startowy

{s1, s2, …, sN} – zbiór N liniowy niezależnych kierunków, np. si = ei



  1. Dokonaj N minimalizacji w kierunku si, używając poprzedniego punktu minimum jako początek kolejnego przeszukiwania. Rozpocznij od kierunku s1 i skończ na sN. Następnie wykonaj jeszcze jedno przeszukiwanie w kierunku s1.

  2. Skonstruuj nowy kierunek sprzężony d korzystając z uogólnionej własności równoległych podprzestrzeni.

  3. Jeżeli ||d|| jest małe lub kierunki poszukiwań liniowo zależne to skończ algorytm, wpp. zastąp sj = sj-1 (dla j = 2, …, N). Zaś s1 = d / ||d|| i  pkt. 1

Uwagi:


- Metoda Powell’a jest powszechnie wykorzystywaną metodą w sytuacji niedostępności informacji o pochodnych

- Metoda ta posiada dowody zbieżności dla funkcji kwadratowych, w przeciwieństwie np. do metody Nelder’a-Mead’a



3.2 Metody gradientowe


Definicja

Kierunek sk nazywamy kierunkiem poprawy (spadku) w punkcie xk, gdy zachodzi nierówność: grad[f(xk)] sk ≤ 0.

Geometryczna interpretacja oznacza, że kierunek poprawy tworzy z kierunkiem minus gradientu kąt mniejszy niż π / 2.

3.2.1 Metoda najszybszego spadku Cauchy’ego3


Dane:

f(x) – funkcja N zmiennych

x0 – punkt startowy

ε – minimalna odległość między kolejnymi wygenerowanymi punktami, tj. ||xk – xk+1||



  1. Wyznaczenie kierunku przeszukiwań:

sk = -grad[f(xk)]

  1. Zdefiniowanie nowego punktu xk+1, jako funkcji α:

xk+1(α) = xk – α sk

  1. Wybranie takiego α, że zminimalizowane jest f[xk+1(α)], a więc:

α = arg minα f[xk+1(α)]

4. Jeżeli ||xkxk+1|| > ε, to k := k + 1 i wrócić do pkt. 1, wpp. zakończyć algorytm i zwrócić xk+1

Uwagi (za Ostatnin(2005)):


  • Szybkość zbieżności metody Najszybszego Spadku jest niedopuszczalnie mała. Wynika to z faktu, że zmiany zmiennych zależą bezpośrednio od wielkości gradientu, który dąży do zera w otoczeniu punktu minimum oraz nieobecny jest mechanizm przyspieszania ruchu do punktu minimum w końcowych iteracjach.

  • Metoda w każdej iteracji poprawia wartość funkcji celu, tj. f(xk+1) < f(xk)

  • Metoda Cauchy’ego pozwala znacząco zmniejszyć wartość funkcji celu przy ruchu z punktów położonych w znacznych odległościach od punktu minimum i dlatego często wykorzystuje się ją przy realizacji metod gradientowych w charakterze początkowej procedury

3.2.2 metoda Newton’a


Dane:

f(x) – funkcja wielu N zmiennych



x0 – punkt startowy

ε – minimalna odległość między kolejnymi wygenerowanymi punktami, tj. ||xkxk+1||



  1. Wyznaczenie kierunku przeszukiwań:

sk = -hesjan[f(xk)]-1 grad[f(xk)]

  1. Zdefiniowanie nowego punktu xk+1, jako funkcji α:

xk+1(α) = xk – α sk

  1. Wybranie takiego α, że zminimalizowane jest f[xk+1(α)], a więc:

α = arg minα f[xk+1(α)]

  1. Jeżeli ||xkxk+1|| > ε, to k := k + 1 i wrócić do pkt. 1, wpp. zakończyć algorytm i zwrócić xk+1

Metoda Newton’a jest bardzo efektywna, gdy RO znajduje się już relatywnie blisko badanego obecnie punktu xk. W otoczeniu RO badane funkcje celu są bardzo dobrze aproksymowane funkcjami kwadratowych, stąd Metodę Newton’a wykorzystuje się zazwyczaj w ostatnich iteracjach procedur poszukiwania RO. Dla punktów odległych od RO, kierunek przeszukiwań Metody Newton’a nie musi być nawet kierunkiem poprawy, tj. kolejne iteracje Metody Newton’a mogą pogorszać wartość funkcji celu. Warunkiem wystarczającym na to aby sk był kierunkiem poprawy jest półdodatniość macierzy: hesjan[f(xk)]-1. Nawet to nie gwarantuje, że spadek wartości funkcji celu będzie wystarczająco duży. W takim przypadku należy zastosować metodę Cauchy’ego. Kolejna metoda łączy zalety metody Cauchy’ego i Newton’a.

3.2.3 Metoda Markwardt’a


Metod’a Markwardt’a łączy zalety metody Cauchy’ego i Newton’a. Dla punktów odległych od RO, metoda przypomina bardziej metodę Cauchy’ego i wraz z zbliżaniem się do RO adaptacyjnie przekształca się w metodę Newton’a.

Dane:


f(x) – funkcja wielu N zmiennych

x0 – punkt startowy

ε – dopuszczalna wielkość modułu gradientu ||grad[f(xk)]||

λ – wartość wyznaczająca początkową rolę Hesjanu w wyznaczaniu kierunku poprawy

M – maksymalna liczba iteracji



  1. Obliczyć ||grad[f(xk)]||

  2. Jeżeli ||grad[f(xk)]|| < ε lub k > M, to zakończ, wpp. Kontynuuj

  3. sk = -{hesjan[f(xk)] + λk I}-1 grad[f(xk)], zaś xk+1(α) = xk + α sk, gdzie: α = argminα f[xk+1(α)], I – macierz jednostkowa

  4. Jeżeli f(xk+1) < f(xk) wówczas λ = ½ λ i  pkt. 1, wpp. λ = 2λ i  pkt. 3

Tym większa wartość λ tym mniejsza rola Hesjanu w wyznaczaniu kierunku poprawy. Stąd przy dużych wartościach λ, metoda Markwardt’a zachowuje się jak metoda Cauchy’ego. Wraz ze zbliżaniem się do RO, wartość λ maleje i kierunek poprawy bardziej naśladuje ten z Metody Newton’a.
      1. Metoda gradientów sprzężonych Fletcher’a-Reeves’a


Dane:

f(x) – funkcja N zmiennych



x0 – punkt startowy

  1. Policz grad[f(x)] i ustal s0 = - grad[f(x0)]

  2. Dokonaj minimalizacji w kierunku s0 z punktu x0. Oznacz RO jako x1. Oblicz grad[f(x1)]

Zainicjuj k = 1.

  1. Oblicz:

sk = - grad[f(xk)] + {|| grad[f(xk)]|| / || grad[f(xk-1)]||} sk-1

Dokonaj minimalizacji w kierunku sk z punktu xk. RO oznacz jako xk.



Jeżeli zaszedł warunek STOPu, zakończ, wpp. k = k +1 i  pkt. 3

3.2.5 Metoda zmiennej metryki




    1. Zadania


  1. [PD] Dla zadanych funkcji pokaż proces zbiegania4 Metody Cauchy i Newton do RO – wyciągnij wnioski jakościowe:

    1. f(x) = 8x12 + 4x1x2 + 5x22

    2. f(x1, x2) = 100(x2-x12)2 + (1 – x1)2 (funkcja Rosenbrock’a)

  2. [PD][Deb] Dla zadanej funkcji f(x1, x2) wykonaj trzy iteracje metody złotego podziału, w celu oszacowania minimalnego punktu na prostej łączącej punkt [ 3,  4]T z [3, 2]T. Ogranicz przeszukiwanie do obszaru pomiędzy tymi punktami.

f(x1, x2) = 2 + (x12 – x2)2 + x22

  1. [Deb] Sprawdź czy dany kierunek poszukiwań s jest kierunkiem poprawy w punkcie x dla funkcji f:

    1. [LAB] f(x1, x2) = 2x12 + x22 – 2x1x2 + 4, s = (1, 1)T, x = (2, 3)T

    2. f(x1, x2) = x14 + x23 – 2 x12x22 + 10x1 / x22, s = (-1, 2)T, x = (0,1)T

    3. [PD] f(x1, x2) = x13 – 2x1x22 + 4x1 + 10, s = (cos ω, sin ω)T, x = (1, 2)T

    4. [PD] f(x1, x2, x3) = (x1 + 2x2 + 3x3)2 – 5(x12 – x3)2, s = (1, 2, -1)T, x = (1, 0, 1)T

    5. f(x1,x2,x3) = x13x2 + 0.2exp(x22 – x3), s = (-10, -20, 5)T, x = (1, 2, 6)T

  2. [LAB][Deb] Wykonaj dwie iteracji metody pełzającego Simplex’u Nelder’a – Mead’a w celu znalezienia nowego simplex’u dla f(x1, x2) = 10 – x1 + x1x2 + x22, x1 = (0, 2)T, x2 = (0, 0)T, x3 = (1, 1)T. Ponadto, załóż: β = 0,5, γ = 2.

  3. [Deb] Wykonaj dwie iteracje poniższych algorytmów dla zadania minimalizacji funkcji f(x1, x2, x3, x4) = (x1 + 2x2 – 1)2 + 5(x3 – x4)2 + (x2 – 3x3)4 + 10(x1 – x4)4

    1. metoda Cauchy’ego x0 = (2, -1, 0, 1)T

  4. [PD][Deb] Rozważ następujące zadanie minimalizacji funkcji:

f(x1, x2, x3) = x12 + 2x22 + 2x32 + 2x1x2 + 2x2x3 – 3x1

Mając kierunek poszukiwań s1 = (1, -1, -1)T oraz punkty startowe x1 = (1, 0, 1)T, x2 = (0, 1, 0)T, znajdź nowy kierunek poszukiwań s2 sprzężony do s1, używając przy tym własności równoległych podprzestrzeni. Pokaż, że wyznaczony kierunek s2 jest sprzężony do s1 w konteście rozpatywanej funkcji f.



  1. [PD] Zaimplementuj jedną z metod:

    1. Nelder’a-Mead’a

    2. Kierunków sprzężonych Powell’a.

    3. Gradientów sprzężonych (0,5pkt)

    4. DFP (0,5pkt)

  2. [PD] Wykonaj trzy iterację metody Markquardt’a dla funkcji Himmelblau’a:

Min f(x) = (x12 + x2 -11)2 + (x1 + x22 – 7)2.

Przyjmij punkt startowy x0 = (0,0) i λ = 100


3.4 Literatura:


http://www.cse.illinois.edu/iem/

Bazaraa M.S., Sherali H.D., Shetty C.M., Nonlinear Programming – Theory and Algorithms, 3rd Edition, Wiley-Interscience, 2006.

Deb K., Optimization for Engineering Design: Algorithms and Examples, PHI Learning Private Limited, New Delhi, 2009.

Kusiak J., Danielewska-Tułecka A., Oprocha P., Optymalizacja. Wybrane metody z przykładami zastosowań., Wydawnictwo Naukowe PWN, Warszawa 2009.

Ostatnin A., Optymalizacja liniowa i nieliniowa, Politechnika Białostocka, Białystok 2005.

Rao, Singiresu S., Engineering Optimization: Theory and Practice, 4th edition, John Wiley & Sons, 2009.

Stachurski A., Wprowadzenie do optymalizacji, Oficyna Wydawnicza PW, Warszawa 2009.

Stadnicki J., Teoria i Praktyka Rozwiązywania zadań optymalizacji z przykładami zastosowań technicznych, WNT, Warszawa 2006.


3.5 Implementacje algorytmów


function [ optimum,wartosc ] = cauchy (f,x0,epsilon)

%ustalenie wartosci epsilon, w przypadku nie podania przez uzytkowanika

%wektor kolumnowy (!!) x0 musi byc podany, aby algorytm zidentyfikowal

%wymiar x'a

if nargin<3

epsilon=1e-4;

endif;

%kierunek poprawy s:



s=-pochodna(f,x0);

%definiowana jest funkcja jednej zmiennej, tj. parametru a stojącegy przy

%kierunku poprawy s

krok=@(a)feval(f,x0+a*s);

%nowy punkt jest generowany w opraciu o wczesniejszy i optymalny krok

x_next=x0+golden(krok,0)*s;

%dalsze kroki wykonywane sa w petli, dopoki odleglosc miedzy kolejnymi

%wygenenowanymi punktami bedzie mniejsza niz epsilon

while modul(x_next-x0)>epsilon

x0=x_next;

s=-pochodna(f,x0);

krok=@(a)feval(f,x0+a*s);

x_next=x0+golden(krok,0)*s;

endwhile;

%ostatecznie algorytm zwraca dwie wartosci: RO i wart. funkcji w RO

optimum=x_next;

wartosc=feval(f,x_next);

endfunction

function [ optimum,wartosc ] = newton (f,x0,epsilon)

%ustalenie wartosci epsilon, w przypadku nie podania przez uzytkowanika

%wektor kolumnowy (!!) x0 musi byc podany, aby algorytm zidentyfikowal

%wymiar x'a

if nargin<3

epsilon=1e-4;

endif;

%kierunek poprawy s:



s=-inv(hesjan(f,x0))*pochodna(f,x0);

%definiowana jest funkcja jednej zmiennej, tj. parametru a stojącegy przy

%kierunku poprawy s

krok=@(a)feval(f,x0+a*s);

%nowy punkt jest generowany w opraciu o wczesniejszy i optymalny krok

x_next=x0+golden(krok,0)*s;

%dalsze kroki wykonywane sa w petli, dopoki odleglosc miedzy kolejnymi

%wygenenowanymi punktami bedzie mniejsza niz epsilon

while modul(x_next-x0)>epsilon

x0=x_next;

s=-inv(hesjan(f,x0))*pochodna(f,x0);;

krok=@(a)feval(f,x0+a*s);

x_next=x0+golden(krok,0)*s;

endwhile;

%ostatecznie algorytm zwraca dwie wartosci: RO i wart. funkcji w RO

optimum=x_next;

wartosc=feval(f,x_next);

endfunction


3.6 Odpowiedzi




1 ang. Direct Search Methods

2 Jest to odbicie punktu xh względem hiperpowierzchni rozpiętej przez zbiór punktów {xi}i≠h

3 ang. Steepest Descent

4 tj. dla RO i f(RO), policzonych metodami analitycznymi, pokaż ||xk – RO|| i [f(xk) – f(RO)], dla k = 0, 1, 2, ...





©snauka.pl 2016
wyślij wiadomość