TECHNIKA OPTYMALIZACJI laboratorium
Ćwiczenie III: Zadanie programowania nieliniowego bez ograniczeń algorytmy optymalizacji lokalnej

1.Cel ćwiczenia

Celem ćwiczenia jest zapoznanie się z jedną z metod rozwiązywania nieliniowego zadania optymalizacji bez ograniczeń na przykładzie algorytmu optymalizacji lokalnej.

Algorytm optymalizacji lokalnej dla zadania nieliniowej optymalizacji bez ograniczeń pozwala na wyliczenie takiego wektora zmiennych decyzyjnych $\hat{x}$, dla którego funkcja celu f(x) osiąga minimum:

$$\operatorname{}{f\left( x \right) = f(\hat{x})}$$

2. Wykonanie ćwiczenia

2.1 Zadanie testowe numer 1

f(x) = (x₁ − 2)² + (x₂ − 2)²

Punkt startowy x₀ = [5,  3] , f(x₀) = 10,0

Punkt optymalny: $\hat{x} = \left\lbrack 2,0\ ;2,0 \right\rbrack,$ $f\left( \hat{x} \right) = 0$

Wizualizacja poziomicy z kolejnymi iteracjami została przedstawiona dla funkcji fminunc (nieużytej w tym przypadku na zajęciach) , ze względu na czytelność wykresu(ilość iteracji):

(przyjęte kryterium stopu TolX = 10^-2,użyta funkcja fminunc)

x = 2.0000 2.0000 fval = 6.6669e-16 iterations = 2

exitflag = 1 algorithm: 'medium-scale: Quasi-Newton line search'

Rys. 1 Poziomica funkcji z naniesionymi kolejnymi iteracjami (funkcja fminunc, TolX = 10-2 )

Analiza wyników uzyskanych na zajęciach:

Wyniki optymalizacji otrzymane za pomocą programu matlab (funkcje: fminsearch/fminunc):

(przyjęte kryterium stopu TolX = 10^-2)

x = 1.9972 2.0023 fval = 1.2933e-005 iterations = 27

exitflag = 1 algorithm: ’Nelder-Mead simplex direct search’

Rys. 2 Wykres funkcji( TolX = 10-2 )

Rys. 3 Zmiany wartości przy kolejnych iteracjach( TolX = 10-2 )

Wyniki po zmianie kryterium (przyjęte kryterium stopu TolX = 10^-8 ):

x = 2.000 2.000 fval = 5.6324e-18 iterations = 73

exitflag = 1 algorithm: ’Nelder-Mead simplex direct search’

Rys. 4 Zmiany wartości przy kolejnych iteracjach( TolX = 10-8 )

Następnie dokonane zostały niezbędne obliczenia analityczne:

Postać naszej funkcji : f(x) = (x₁ − 2)² + (x₂ − 2)² 

Gradient funkcji:

Wyznaczanie punktu stacjonarnego:

grad f(x)=0 (zależność na punkt stacjonarny)

2x₁ − 4 = 0 →  x₁ = 2

2x₂ − 4 = 0 →  x₂ = 2

Możliwe ekstremum w punkcie: P=(2,2)

Przechodzimy do macierzy Hessego :

Wyznacznik z powyższej macierzy:

Jako , że wyznacznik >0 (4>0) w punkcie (2,2) znajduje się ekstremum lokalne.

Następnie z odpowiedniej zależności:

wnioskujemy , że P(2,2) jest minimum ,gdyż 2>0.

f(2,2) = (2 − 2)² + = (2 − 2)² = 0

Wartość funkcji w minimum to 0(zgodnie z podanymi założeniami).

2.1 Zadanie podane przez prowadzącego

f(x) = (x₁ − 2)² + (x₁ − x₂²)²

Z wykresu wywnioskować można , że funkcja posiada dwa minima leżące do siebie równolegle względem osi y. Do obliczenia gradientów oraz hesjanów wykorzystano funkcje fminunc . Oba obliczenia zostały wykonane dla dwóch minimów przy TolX = 10^-8)

I minimum

Punkt startowy x₀ = [2  1.6]

Wyniki optymalizacji otrzymane za pomocą programu matlab (funkcje: fminsearch/fminunc):

(przyjęte kryterium stopu TolX = 10^-2)

Wizualizacja poziomicy z kolejnymi iteracjami została przedstawiona dla funkcji fminunc również ze względu na większą czytelność wykresu(ilość iteracji):

x = 2.0476 1.4313 fval = 0.0023 iterations = 2

exitflag = 1 algorithm: 'medium-scale: Quasi-Newton line search'

Rys.5 Poziomica funkcji z naniesionymi kolejnymi iteracjami ( funkcja fminunc,TolX = 10-2)

Analiza wyników uzyskanych na zajęciach:

x = 2.0033 1.4133 fval = 4,5044e-005 iterations = 19

exitflag = 1 algorithm: ’Nelder-Mead simplex direct search’

Rys. 6 Wykres funkcji( TolX = 10-2 )

Rys. 7 Zmiany wartości przy kolejnych iteracjach( TolX = 10-2 )

Wyniki po zmianie kryterium (przyjęte kryterium stopu TolX = 10^-8 ):

x = 2.0000 1.4142 fval = 3,6486e-017 iterations = 65

exitflag = 1 algorithm: ’Nelder-Mead simplex direct search’

Wartości zaczerpnięte z zadania(funkcja fminunc):

Otrzymany gradient: $\mathbf{\nabla}\mathbf{f}\mathbf{(}\hat{\mathbf{x}\mathbf{)}}\mathbf{=}$ 1.0e-06*$\begin{bmatrix} \mathbf{\text{\ \ }}\mathbf{0}\mathbf{.}\mathbf{1381} \\ \mathbf{\ -}\mathbf{0}\mathbf{.}\mathbf{3026} \\ \end{bmatrix}$

$\nabla f(\hat{x)}\ \cong 0$ : możemy wywnioskować , że spełniony jest warunek stacjonarności

Otrzymany hesjan: H =$\begin{bmatrix} \mathbf{\ }\mathbf{4}\mathbf{.}\mathbf{000} & \mathbf{- \ }\mathbf{5}\mathbf{.}\mathbf{6572} \\ \mathbf{\ \ -}\mathbf{5}\mathbf{.}\mathbf{6572} & \mathbf{\ 16.0059} \\ \end{bmatrix}$

Hesjan jest dodatnio określony:
det H ≅ 32,0196 >0
oraz zauważamy , że 4.000>0

W naszym punkcie znajduje się ekstremum , jest ono minimum.

Rys. 8 Zmiany wartości przy kolejnych iteracjach( TolX = 10-8 )

II minimum

Punkt startowy x₀ = [2    − 1.6]

Wyniki optymalizacji otrzymane za pomocą programu matlab (funkcje: fminsearch/fminunc):

(przyjęte kryterium stopu TolX = 10^-2)

x = 2.0033 -1.4133 fval = 4,5044e-005 iterations = 19

exitflag = 1 algorithm: ’Nelder-Mead simplex direct search’

Rys. 9 Zmiany wartości przy kolejnych iteracjach( TolX = 10-2 )

Zmiana kryterium stopu TolX = 10^-8:

x = 2.0000 -1.4142 fval = 3,6486e-017 iterations = 65

exitflag = 1 algorithm: ’Nelder-Mead simplex direct search’

Wartości zaczerpnięte z zadania(funkcja fminunc):

Otrzymany gradient: $\mathbf{\nabla}\mathbf{f}\mathbf{(}\hat{\mathbf{x}\mathbf{)}}\mathbf{=}$ 1.0e-06*$\begin{bmatrix} \mathbf{\text{\ \ }}\mathbf{0}\mathbf{.}\mathbf{1381} \\ \mathbf{\text{\ \ }}\mathbf{0}\mathbf{.}\mathbf{3026} \\ \end{bmatrix}$

$\nabla f(\hat{x)}\ \cong 0$ : możemy wywnioskować , że spełniony jest warunek stacjonarności

Otrzymany hesjan: H =$\begin{bmatrix} \mathbf{\ }\mathbf{4}\mathbf{.}\mathbf{000} & \mathbf{\ }\mathbf{5}\mathbf{.}\mathbf{6572} \\ \mathbf{\text{\ \ \ }}\mathbf{5}\mathbf{.}\mathbf{6572} & \mathbf{\ 16.0059} \\ \end{bmatrix}$

hesjan jest dodatnio określony, ponieważ wszystkie wiodące minory główne są dodatnie::
det H ≅ 32,0196 >0

oraz zauważamy , że 4.000>0

W naszym punkcie znajduje się ekstremum , jest ono minimum.

Rys. 10 Zmiany wartości przy kolejnych iteracjach( TolX = 10-8 )

3.Opis algorytmu

Metoda pełzającego simpleksu Neldera-Meada

Metoda MelderaMetoda ta polega na utworzeniu w przestrzeni Eⁿ⁺¹ n-wymiarowego simplexu o n+1 wierzchołkach tak, aby można było go wpisać w powierzchnię reprezentującą badaną funkcję celu. Jednowymiarowym simplexem jest odcinek o dwóch wierzchołkach, simplexem dwuwymiarowym jest trójkąt i ogólnie simplexem n-wymiarowym o n+1 wierzchołkach jest zbiór wszystkich punktów określonych przez wektory:

czyli jest to wielościan o n+1 wierzchołkach rozpiętych na n+1 wektorach bazowych (S_j). Współrzędne punktów simplexu oznaczono jako x_j. 
Na początku procedury wylicza się współrzędne punktów wierzchołkowych simplexu P_j (dla j = 1 .. n+1) przy założeniu pewnej odległości między tymi wierzchołkami (czyli kroku). W następnych iteracjach dokonuje się przekształceń simplexu aż odległość pomiędzy jego wierzchołkami w pobliżu poszukiwanego ekstremum będzie mniejsza od założonej dokładności obliczeń e. To właśnie zostało przyjęte jako kryterium zbieżności dla tej metody. 
Podczas wykonywania algorytmu metody simplex stosuje się trzy podstawowe operacje:

- odbicie punktu P_h względem P': P* = (1 + a)P' - aP_h

- ekspansja punktu P** względem P': P** = (1 - c)P* - cP'

-kontrakcja punktu P_h względem P': P*** = bP_h + (1 - b)P'

Oznaczenia:

x₀ - punkt startowy 
e - wymagana dokładność obliczeń 
P_h - wybrany punkt wierzchołkowy simplexu spośród n+1 wierzchołków P_i, w którym wartość badanej funkcji osiąga maksimum. 
P_L - wybrany punkt wierzchołkowy simplexu spośród n+1 wierzchołków P_i, w którym wartość badanej funkcji osiąga minimum. 
P' - środek symetrii simplexu z wyłączeniem punktu P_h zdefiniowany jako:

d - początkowa odległość pomiędzy wierzchołkami wyjściowego simplexu 
a - współczynnik odbicia (a>0) 
b - współczynnik kontrakcji (0<b<1) 
c - współczynnik ekspansji (c>1) 
n - liczba zmiennych niezależnych

Wartości współczynników a, b, c dobiera się w sposób eksperymentalny, choć w przykładach rozpatrywanych przez autorów metody jako optymalne przyjęto wartości a = 1, b = 0,5 oraz c = 2.

Algorytm:

1) Obliczenie wartości funkcji celu w punktach wierzchołkowych simplexu F_j = f(P_j) dla j = 1 .. n+1 
2) Wyznaczenie h i L takich, że f(P_h) = max oraz f(P_L) = min spośród zbioru F_j. 
3) Obliczenie środka symetrii simplexu P'. 
4) Odbicie P* punktu P_h względem P'. 
5) Obliczenie wartości funkcji F_s = f(P') oraz F_o = f(P*).

Jeśli F_o < min, to: 
6) Obliczenie P** (ekspansja) i wartości funkcji F_e = f(P**). 
7) Jeśli F_e < max, to podstawiamy P_h = P**, w przeciwnym przypadku podstawiamy P_h = P*. 
8) Powtórzenie procedury od kroku 1), o ile nie jest spełnione kryterium na minimum.

Jeśli F_o > min, to: 
6) Jeśli F_o >= f(P_j) dla j = 1 .. n+1 (z pominięciem j = h) i F_o >= max, przejście do następnego kroku. 
    Jeśli zaś F_o < max, to podstawiamy P_h = P*. 
7) Wykonanie kontrakcji P*** punktu P_h względem P'. 
8) Obliczenie F_k = f(P***). 
9) Jeżeli F_k >= max to wykonujemy redukcję simplexu wg wzoru: P_j = 0,5*(P_j + P_L) dla j = 1 .. n+1. 
    Jeżeli natomiast F_k < max, to podstawiamy P_h = P*** i przechodzimy do kroku 11). 
10) Jeśli F_o < f(P_j) dla j = 1 .. n+1 (z pominięciem j = h), podstawiamy P_h = P*. 
11) Powtórzenie procedury od kroku 1), jeśli nie zostało spełnione kryterium na minimum. 

Źródło: http://www.kmg.zut.edu.pl/opt/wyklad/bezgrad/simplex.html

Metoda Quasi-Newtona (zmiennej metryki)

Metoda Quasi-Newtona może być używana, gdy obliczenie Hessianu (macierzy drugich pochodnych funkcji) wymaganego przez metodę Newtona  jest trudne lub czasochłonne. Idea metody polega na przybliżeniu Hessianu lub jego odwrotności za pomocą pierwszych pochodnych.

W każdym kroku algorytmu przybliżenie Hessianu (H) lub jego odwrotności (B) jest uaktualniane przy pomocy informacji o gradiencie. Informacje o krzywiźnie funkcji są gromadzone w trakcie wykonywania obliczeń.

Oznaczenia:

H – Hessian

B – odwrotność macierzy H – 

x₀ – pierwsze przybliżenie rozwiązania (punkt startowy)

x_i – i-te przybliżenie rozwiązania

H₀, B₀ – macierze początkowe

H_i, B_i – i-te macierze przybliżone

H_i^u, B_i^u – i-te poprawki wartości macierzy H_i, B_i

ε – wymagana dokładność

Metody przybliżania wartości H i B:

Istnieją różne metody obliczania przybliżeń, z których dwie najpopularniejsze to:

- DFP (Davidon-Fletcher-Powell)

-BFGS (Broyden-Fletcher-Goldfarb-Shanno)

co po odwróceniu daje:

Ponadto istnieje cała rodzina metod znana jako Broyden family:

, gdzie Φ jest dowolną liczbą rzeczywistą

Algorytm:

1. Ustal i = 0, x₀, B₀ oraz kryterium zatrzymania ε.

, gdzie macierz początkowa H₀ może być dowolną symetryczną, dodatnie określoną macierzą. W szczególności .

2. Sprawdź, czy punkt x_i spełnia kryterium stopu – np. jeślito x_i jest rozwiązaniem

3. Wyznacz kierunek poprawy

4. Wyznacz kolejne przybliżenie rozwiązania , gdzieto długość kroku minimalizująca jednowymiarową funkcję

5. Wyznacz poprawkęzgodnie z wybraną metodą (DFP, BFGS), używając wartości , oraz 

6. Uaktualnij

7. i = i+1

8. Idź do punkt 2.

Źródło: http://optymalizacja.w8.pl/QuasiNewton.html

4.Wnioski

Metoda z której korzysta funkcja „fminunc” wykonuje dużo mniej iteracji niż metoda Neldera-Meada używana przez funkcję „fminsearch”.Znaczący zysk w ilości iteracji jest jednak okupiony kosztami. Dodatkowym kosztem użycia metody Quasi-Newtonoskiej jest zwiększona potrzeba obliczeniowa. W każdej z iteracji obliczany jest hessian oraz gradient – czego nie zauważamy przy Nelderze-Meadzie. Metoda Quasi-Newtonoska pozwala na aproksymację hesjanu za pomocą pierwszych pochodnych , co ułatwia obliczenia z nim związane. Eksperymenty związane ze zmianą kryterium stopu TolX z 10^-8 na 10^-2 pozwalają na wyciagnięcie dwóch wniosków. Zmiana tego kryterium wpływa na ilość iteracji (znacząco mniejsza ilość przy kryterium ustawionym na 10^-2 ) . Również wynik jest przedstawiony z różną dokładnością , co ma związek z ilością iteracji oraz długością obliczeń wykonywanych przez nasz algorytm.