TECHNIKA OPTYMALIZACJI - LABOLATORIUM
Ćwiczenie III: Zadanie programowania nieliniowego bez ograniczeń
algorytmy optymalizacji lokalnej.
Skład grupy:
Prowadzący:
Data wykonania
ćwiczenia:
- Dawid Cekus
- Mateusz Grzelak
Dr inż. Ewa
Szlachcic
3.12.2014
1. Cel ćwiczenia
Celem ćwiczenia było poznanie metody rozwiązywania nieliniowego zadania
optymalizacji. Na funkcję celu nie zostały nałożone żadne ograniczenia.
Wykorzystany został algorytm optymalizacji lokalnej. Algorytm ten umożliwia
wyznaczenie takiego wektora zmiennych decyzyjnych x, dla którego funkcja celu
f(x) osiąga minimum:
min f (x )=f ( ^x)
x ∈R
n
2. Przebieg ćwiczenia
2.1. Zadanie testowe nr 6
Funkcja dla której będzie poszukiwana wartość minimalna ma postać:
f (x
1
, x
2
)=
x
1
4
+
x
2
4
−
2x
1
2
x
2
−
4 x
1
+
3
W celu obliczenia minimum lokalnego posłużono się narzędziem fminunc.
Wyniki działania procedury dla wartości TolFun=0.01 oraz punktu startowego
x0=[1.5;1]:
x=[1.3129;0.9542]
fval=-1.7410
exitflag=1
output:
iterations:3
algorithm: 'medium-scale: Quasi-Newton line search'
hessian=
(
16.8724 −5.2519
−
5.2519 10.9298
)
grad=
(
0.0408
0.0283
)
Wykres funkcji celu:
Wykres warstwicy funkcji wraz z punktami uzyskanymi w kolejnych iteracjach dla
TolFun=0.01:
Wykres przedstawiający wartości funkcji celu w kolejnych iteracjach dla
TolFun=0.01
Wyniki działania procedury dla wartości TolFun=0.00001 oraz punktu startowego
x
0
=[1.5;1]:
x=[1.3091;0.9498]
fval=-1.7411
exitflag=1
output:
iterations:7
algorithm: 'medium-scale: Quasi-Newton line search'
hessian=
(
16.7702 −5.2367
−
5.2367 10.8284
)
grad=
1.0e-5∗
(
−
0.1184
−
0.1431
)
Wykres warstwicy funkcji wraz z punktami uzyskanymi w kolejnych iteracjach dla
TolFun=0.00001:
Wykres przedstawiający wartości funkcji celu w kolejnych iteracjach dla
TolFun=0.00001
Rozwiązanie analityczne:
f (x
1
, x
2
)=
x
1
4
+
x
2
4
−
2x
1
2
x
2
−
4 x
1
+
3
∇
f (x
1,
x
2
)=[
4x
1
3
−
4x
1
x
2
−
4,4 x
2
3
−
2x
1
2
]
Warunek konieczny istnienia ekstremum w punkcie mówi, że gradient w tym punkcie
jest równy zero:
∇
f (x
1,
x
2
)=
0
W przypadku funkcji numer 6 warunek ten jest spełniony dla:
x
1
=
1.309083347
x
2
=
0.9498059087
Wyznaczony punkt jest określany jako punkt stacjonarny. Może się tam znajdować
ekstremum. W celu zbadania czy w powyższym punkcie znajduje się ekstremum
należy zbudować macierz drugich pochodnych tzw. macierz Hessego(hesjan).
W tym przypadku hesjan ma nastepującą postać ogolną:
(
12x
1
2
−
4 x
2
−
4 x
1
−
4 x
1
12 x
2
2
)
W punkcie stacjonarnym ma następującą postać:
(
16.765166877912669 −5.236333388000000
−
5.236333388000000 10.825575170417192
)
Macierz jest dodatnio określona. Świadczy to o tym, że w punkcie stacjonarnym
znajduje się minimum lokalne.
2.2. Zadanie testowe nr 17
Funkcja dla której będą poszukiwane minima lokalne oraz globalne ma postać:
f (x
1,
x
2
)=
4 x
1
2
−
2.1 x
1
4
+
1
3
x
1
6
+
x
1
x
2
−
4 x
2
2
+
4 x
2
4
Wykres funkcji celu:
Wykres warstwicowy funkcji celu:
Na podstawie tego wykresu widać, że funkcja posiada 6 minimów(4 lokalne i 2
globalne) i pomaga dobrać punkty startowe dla poszukiwań. Dodatkowo minima są
pogrupowane w 3 pary. Minima z każdej z par są symetryczne względem punktu
(0,0).
W celu odnalezienia wszystkich minimów posłużono się narzędziem fminsearch.
Pierwsze minimum obliczone dla punktu startowego x
0
=[0,1] dla TolFun=0.1
x
0
=[0;1]:
x=[-0.0898;0.7127]
fval=-1.0316
exitflag=1
output:
iterations:48
algorithm: 'Nelder-Mead simplex direct search'
Procedura fminsearch nie posiada możliwości policzenia gradientu ani hesjanu w
punkcie dlatego obliczenia te zostały przeprowadzone przy użyciu narzędzi
dostępnych w programie Maple.
∇
f (−0.0898,0 .7127)=
(
0.0003711875
0.00075611813
)
H (−0.0898,0 .7127)=
(
7.797436479
1
1
16.38118192
)
Wykres warstwicy funkcji wraz z punktami uzyskanymi w kolejnych iteracjach dla
TolFun=0.1 (pierwsze minimum):
Wykres przedstawiający wartości funkcji celu w kolejnych iteracjach dla
TolFun=0.1(pierwsze minimum):
Drugie minimum zostało odnalezione dla punktu startowego x
0
=[-1.75,0.8] oraz dla
TolFun=0.1
x
0
=[-1.75;0.8]:
x=[-1.7036;0.7961]
fval=-0.2155
exitflag=1
output:
iterations:23
algorithm: 'Nelder-Mead simplex direct search'
∇
f (−1.7036, 0.7961)=
(
0.00014466
0.000375115
)
H (−0.0898,0 .7127)=
(
19.09394785
1
1
22.42121008
)
Wykres warstwicy funkcji wraz z punktami uzyskanymi w kolejnych iteracjach dla
TolFun=0.1 (drugie minimum):
Wykres przedstawiający wartości funkcji celu w kolejnych iteracjach dla
TolFun=0.1(drugie minimum):
Trzecie minimum zostało odnalezione dla punktu startowego x
0
=[-.175,-0.6] oraz dla
TolFun=0.1
x
0
=[-1.75;0.8]:
x=[-1.6071;-0.5686]
fval=2.1043
exitflag=1
output:
iterations: 25
algorithm: 'Nelder-Mead simplex direct search'
∇
f (−1.6071,−05686)=
(
0.00009718
0.000391698
)
H(−1.6071,−0.5686)=
(
9.62121558
1
1
7.51868608
)
Wykres przedstawiający wartości funkcji celu w kolejnych iteracjach dla
TolFun=0.1(trzecie minimum):
Czwarte minimum zostało odnalezione dla punktu startowego x
0
=[0,-1] oraz dla
TolFun=0.0001
x
0
=[0;-1]:
x=[0.0898;-0.7127]
fval= -1.0316
exitflag=1
output:
iterations: 49
algorithm: 'Nelder-Mead simplex direct search'
∇
f (0.0898,−0.7127)=
(
−
0.0003711875
−
0.00075611813
)
H(0.0898,−0.7127)=
(
7.797436479
1
1
16.38118192
)
Wykres przedstawiający wartości funkcji celu w kolejnych iteracjach dla
TolFun=0.0001(czwarte minimum):
Piąte minimum zostało odnalezione dla punktu startowego x
0
=[1.75,-0.8] oraz dla
TolFun=0.0001
x
0
=[1.75;-0.8]:
x=[1.7036;-0.7961]
fval=
exitflag=1
output:
iterations: 23
algorithm: 'Nelder-Mead simplex direct search'
∇
f (1.7036,−0.7961)=
(
−
0.00014466
−
0.000375115
)
H(1.7036,−0.8)=
(
19.09394785
1
1
22.42121008
)
Szóste minimum zostało odnalezione dla punktu startowego x
0
=[1.75,-0.6] oraz dla
TolFun=0.0001
x
0
=[1.75;0.6]:
x=[1.6071,0.5686]
fval= 2.1043
exitflag=1
output:
iterations: 25
algorithm: 'Nelder-Mead simplex direct search'
∇
f (1.6071,0 .5686)=
(
−
0.000009718
−
0.0000391698
)
H(1.6071,0.5686)=
(
9.62121558
1
1
7.51868608
)
Zarówno dla piątego minimum jak i szóstego wykres wartości funkcji celu w
zależności od iteracji był identyczny jak dla minimów symetrycznych względem
punktu (0,0)
3. Opis wykorzystanych metod optymalizacji.
Podczas zajęć były wykorzystywane dwa narzędzia pakietu Matlab:
fminsearch oraz fmiunc. Pierwsze z nich wykorzystuje metodę pełzającego simplexu
Neldera-Meada natomiast drugie korzysta z metod Quasi-newtonowskie.
3.1 Metoda pełzającego simplexu Neldera-Meada.
Metoda pełzającego simplexu jest metodą bezgradientową. W tej metodzie
przeszukuje się kolejno pewne dane kierunki przy wykorzystaniu wersorów osi
współrzędnych.
W tej metodzie działanie jest rozpoczynane od uporządkowania (n+1) punktów
x
1
,x
2
,...,x
n+1
w następujący sposób:
f
n+1
⩾
f
n
⩾
...⩾f
2
⩾
f
1
Wierzchołki te tworzą sympleks o (n+1) wierzchołkach.
W każdej iteracji punkt x
n+1
jest zastępowany punktem „lepszym”
Punkt ciężkości c jest zdefiniowano następująco:
c=
1
n
∑
i=1
n
x
i
Przebieg jednej iteracji metody ma następującą postać:
Konstruowany jest punkt próbny, który jest odpowiedzialny za odbicie punktu
x
n+1
względem hiperpowierzchni rozpiętej na punktach x
1
, .., x
2
, który ma postać:
x
r
=
c+α(c−x
n+1
)
,α>0
α jest współczynnikiem odbicia. Wartość funkcji w punkcie x
r
jest oznaczana: f
r
.
Operacja konstrukcji punktu próbnego jest tożsama z przesunięciem z punktu c w
kierunku (c - x
n+1
) ze współczynnikiem α.
Należy dokonać analizy położenia wartości f
r
.
a) f
1
<= f
r
<= f
n
Punkt x
n+1
jest zastępowany przez x
r
i następuje przejście do następnej iteracji.
Odrzucany jest najgorszy z poprzednich punktów. Może istnieć konieczność
ponownego uporządkowania ze względu na wartość funkcji.
b) fr < f1
Uzyskany punkt jest lepszy od najlepszego dotychczasowych.
Dokonywana jest próba ekspansji(rozciągnięcia) w tym kierunku. Konstruowany jest
następujący punkt:
x
e
=
x+β(x
r
−
c),β>1
β- współczynnik ekspansji
Jeżeli f
e
< f
r
to x
n+1
=x
e
w przeciwnym wypadku x
n+1
=x
r
. Punkty należy znów
uporządkować tak aby wartości funkcji były uporządkowane od najmniejszej do
największej.
c)f
r
> f
n
W tym przypadku wartość próbna może być porównywalna tylko z najgorszą.
W zaistniałej sytuacji istnieje wymóg aby zmniejszyć simpleks. Spowoduje to fakt,
iż będzie można liczyć na ewentualną poprawę wartości funkcji. Należy wówczas
wykonać operację kontrakcji. Należy dokonać wyboru dwóch sympleksów leżących
albo po jednej albo po drugiej stronie hiperpowierzchni, która jest rozpięta na
punktach x
1
, …, x
n
.
Kryterium wyboru odnosi się w tym przypadku do wzajemnej relacji wartości funkcji
w wierzchołku (n+1) -f
r
i f
n+1
gdy f
r
>= f
n+1
x
e
=
c+γ(x
n+1
−
c)
gdy f
r
< f
n+1
x
e
=
c+γ(x
r
−
c)
γ- współczynnik kontrakcji
0 < γ < 1
Jeżeli f
e
< min{f
r
, f
n+1
}, wtedy x
n+1
= x
e
następuje wówczas przejście do kolejnego
kroku. W przypadku gdy ten warunek jest niespełniony operacja ściśnięcia.
Operacje ściśnięcia można zastąpić ściśnięciem całego sympleksu wokół punktu
dotychczas najlepszego
x
i
=
x
1
+γ(
x
i
−
x
1
)
,i=2,...,n+1
Należy również pamiętać o przywróceniu uporządkowania punktów ze względu na
wartość funkcji
Metoda sympleksu Neldera-Meada jest stosowana w odniesieniu do funkcji, które
mogą być nieciągłe albo nieróżniczkowalne. Skuteczność metody ogranicza się do
zadań o małym wymiarze.
Liczba iteracji konieczna do odnalezienia minimum jest większa niż w przypadku
zastosowania metod wykorzystujących pochodne.
3.2 Metody Quasi-newtonowskie
Metody quasi-newtonowskie stanowią alternatywę wobec dyskretnej metody
Newtona. Metoda wykorzystuje informacje o pochodnych drugiego rzędu bez
konieczności jawnego wyliczania macierzy hesjanu. Wykorzystywana jest informacja
o pochodnych pierwszego rzędu w aktualnym punkcie x
k
oraz punkcie następnym
x
k+1
. Wyrażenia opisujące pochodne pierwszego i drugiego rzędu na kierunku s
k
są
uzyskiwane na podstawie rozwinięcia gradientu w szereg Taylora względem
przyrostu zmiennych niezależnych. Wspomniane wyżej wyrażenia mają następującą
postać:
~f'(0)=∇ f (x
k
)
T
s
k
=
g(x
k
)
T
s
k
~f' '(0)=s
kT
∇
2
f (x
k
)
s
k
=
s
kT
G
k
s
k
Wyrażenie, które pozwala określić krzywiznę minimalizowanej funkcji f wzdłuż
kierunku s
k
:
s
kT
G
k
s
k
≈(
g(x
k
+
s
k
)−
g(x
k
))
T
s
k
Metody quasi-newtonowskie wykorzystują fakt, że skoro można oszacować
krzywiznę wzdłuż kierunku na podstawie różnic pochodnych pierwszego rzędu to
można dokonywać próby aproksymacji hesjanu na podstawie informacji
analogicznej. W metodach tego typu aproksymujemy macierz G
k
bądź też jej
odwrotność (G
k
)
-1
. Na początku k-tej iteracji do dyspozycji B
k
, H
k
. W pierwszym
przypadku kolejny kierunek poszukiwań jest wyznaczony z następującej zależności:
B
k
d
k
=−
g
k
w drugim przypadku :
d
k
=−
H
k
g
k
W przypadku braku informacji przyjmuje się H
0
=I(B
0
=I).
I-macierz identycznościowa
Po znalezieniu kolejnego przybliżenia rozwiązania- punktu x
k+1
macierz H
k
lub B
k
jest
modyfikowana w następujący sposób:
H
k+1
=
H
k
+
U
k
gdzie U
k
jest poprawką wyznaczoną na podstawie przyrostów zmiennych
niezależnych i przyrostów gradientu. Poprawka jest konstruowana jest w taki sposób,
aby kolejna macierz spełniała warunek quasi-newtonowski:
H
k+1
r
k
=
s
k
W metodach quasi-newtonowskich rozważa się w głównie poprawki pierwszego i
drugiego rzędu.
Opis metod został stworzony na podstawie książki:
„Podstawy optymalizacji” A.Stachurski, A.P Wierzbicki.
4. Wnioski
Dla pierwszego zadania testowego zostało znalezione jedno minimum lokalne. W
celu jego znalezienia skorzystano z narzędzia fminunc,które korzysta z metod
quasi-newtonowskich. Obliczenia dokonano dla dwóch wartości parametru
TolFun=0.01 oraz TolFun=0.00001. Dla mniejszej wartości parametru można
dostrzec, że wartość rozwiązania jest bardziej dokładna, dodatkowo można
zaobserwować większą liczbę iteracji.
Dla drugiego zadania testowego znaleziono 6 minimów- 2 globalne, 4 lokalne.
Wykreślenie warstwic funkcji było pomocne dla określenia poszczególnych punktów
startowych od których rozpoczną się poszukiwania. W tym przypadku skorzystano z
narzędzia fminsearch, które korzysta z metody pełzającego sympleksu Neldera-
Meada. Pierwsze 3 minima zostały znalezione dla wartości parametru TolFun=0.1 a
następne 3 dla TolFun=0.0001. Zmiana wartości TolFun nie wpłynęła na poprawę
dokładności wyniku oraz na liczbę iteracji.
Znalezione minima są pogrupowane w trzy pary. Minima z każdej pary są
symetryczne względem punktu (0,0).
Na podstawie wyników obliczeń można dostrzec, że sympleks Neldera-Meada jest
mniej efektywny niż metody quasi-newtonowskie. Świadczy o tym liczba iteracji
potrzebna do znalezienia szukanego punktu.