ˇ sen´ı u´kolu z 9. pˇredn´aˇsky (DU10) Reˇ X35ORR Petr Svoboda,
[email protected] 7. ledna 2009
Pˇ r´ıklad 1 √
13 na 32 desetinn´ ych m´ıst je √ 13 = 3, 60555127546398929311922126747049.
Mnohem zaj´ımavˇejˇs´ı, neˇz samotn´ y v´ ysledek, je zp˚ usob v´ ypoˇctu. Pro hledan´e x = plat´ı g(x) = x2 − 13.
(1) √
13 na 32 (2)
Z obecn´eho vzorce pro Newtonovu metodu xk+1 = xk −
g(xk ) g 0 (xk )
(3)
odvod´ım konkr´etn´ı iteraˇcn´ı vzorec xk+1
x2k − 13 = xk − . 2xk
(4)
ˇ Tento vzorec velice rychle konverguje1 ke spr´avn´emu v´ ysledku. Podle slov profesora Stechy, s kaˇzdou dalˇs´ı iterac´ı se zdvojn´asobuje poˇcet platn´ ych m´ıst v´ ysledku. Pochopitelnˇe ale z´aleˇz´ı na poˇca´teˇcn´ı podm´ınce iteraˇcn´ı rovnice x0 . Pro nˇekolik hodnot poˇca´teˇcn´ı podm´ınky x0 jsem vynesl poˇcet platn´ ych m´ıst v z´avislosti na poˇctu iterac´ı rovnice (4). V´ ysledky jsou vidˇet v tabulce 1. Je vidˇet, ˇze zhruba plat´ı to, ˇze se s kaˇzdou iterac´ı zdvojn´asobuje poˇcet desetinn´ ych m´ıst. Pr˚ ubˇeh vzd´alenosti odhadu optima xk od skuteˇcn´eho optima x∗ je velice z´avisl´ y na poˇca´teˇcn´ı podm´ınce. Z toho plyne, ˇze pokud by jsme pˇredem napevno“ stanovili, ” kolik provedeme iterac´ı rovnice (4), tak pro vhodnou volbu poˇca´teˇcn´ı podm´ınky by jsme iterovali zbyteˇcnˇe dlouho. Naopak pro nevhodnou volbu poˇc´ateˇcn´ı podm´ınky by jsme pro pˇredem urˇcen´ y poˇcet iterac´ı nemuseli dos´ahnout poˇzadovan´e pˇresnosti. Program ˇreˇs´ıc´ı tento pˇr´ıklad je v pˇr´ıloze A. Program je napsan´ y v jazyce PHP, protoˇze jsem potˇreboval nˇejak´ y n´astroj s obecnˇe ˇra´dovou aritmetikou. PHP toto umoˇzn ˇuje tak, ˇze ˇc´ısla jsou reprezentov´ana textov´ ymi ˇretˇezci a k operac´ım s nimi jsou pouˇzity funkce bcadd, bcmull, bcdiv atd . . . 1
. . . s ˇr´ adem konvergence 2.
1
poˇcet platn´ ych m´ıst iterace x0 = 1 x0 = 4 x0 = 100 1 0 3 0 2 0 5 0 3 3 11 0 4 4 21 0 5 8 41 0 6 15 83 3 7 32 165 4 8 64 329 8 9 127 659 17 10 253 1317 32 11 506 2631 65 √ Tabulka 1: Poˇcet desitinn´ ych m´ıst ˇc´ısla 13 z´ıskan´eho Newtonovou iteraˇcn´ı metodou.
Pˇ r´ıklad 2 (a) Nalezen´ı min{cos(x)}, x ∈ (0, 2π)
(5)
x
pomoc´ı Newtonovy metody spoˇc´ıv´a v tom, ˇze si funkci cos(x) v bodˇe xk aproximuji kvadratickou funkc´ı g(x)|xk = cos(xk ) − sin(xk )(x − xk ) −
1 cos(x − xk )2 . 2
(6)
Bod xk je m˚ uj aktu´aln´ı odhad optima. Iterativnˇe jej posouv´am podle xk+1 = xk −
g(xk ) dg(xk ) dxk
.
(7)
Kritick´ y krok je pˇritom aproximace cos(x) funkc´ı g(x). Situace je zn´azornˇena na obr´azku 1. 1.5
1.5 cos(x) g(x) x* xk xk+1
1
1.5 cos(x) g(x) x* xk xk+1
1
0.5
0.5
0.5
0
0
0
−0.5
−0.5
−0.5
−1
−1
−1
−1.5
0
1
2
3
4
5
6
−1.5
0
1
2
3
x
4 x
5
6
cos(x) g(x) x* xk xk+1
1
−1.5
0
1
2
3
4
5
6
x
Obr´azek 1: Aproximace cos(x) kvadratickou funkc´ı v r˚ uzn´ ych bodech. Pro body xk bl´ızk´e minimu x∗ se odhad minima xk+1 pomoc´ı kvadratick´a funkce g(x) pˇr´ıliˇs neliˇs´ı od skuteˇcn´eho minima. Se zvˇetˇsuj´ıc´ı vzd´alenost´ı xk od x∗ se kvalita aproximace pomoc´ı 2
, 2π) se parabola g(x) u ´plnˇe otoˇc´ı a vede algoritg(x) zhorˇsuje, aˇz pro hodnoty xk ∈ (0, π2 i ∪ h 3π 2 mus pˇresnˇe opaˇcn´ ym smˇerem. V t´eto mnoˇzinˇe hodnot je pomysln´a Hassova matice2 negativnˇe 2 f (x) semidefinitn´ı, tedy druh´a derivace ∂ ∂x ≤ 0. Nen´ı tedy pˇrekvapiv´e, ˇze Newtonova metoda nefunguje na cel´em intervalu (0, 2π), protoˇze funkce cos(x) nen´ı na cel´e t´eto mnoˇzivˇe konvexn´ı. Zaj´ımav´e je ale to, ˇze pokud by jsme se omezili na obor x ∈ ( π2 , 3π ), na kter´em je cos(x) 2 konvexn´ı, nebude st´ale Newtonova metoda konvergovat pro vˇsechna poˇc´ateˇcn´ı x. Situace je dobˇre vidˇet na prostˇredn´ım obr´azku 1. Zde je evidentnˇe xk > π2 a pˇritom n´asleduj´ıc´ı krok xk+1 je vzd´alenˇejˇs´ı od minima x∗ neˇz poˇc´ateˇcn´ı bod xk . Simulacemi jsem urˇcil, ˇze algoritmus konverguje v oblasti poˇca´teˇcn´ıch podm´ınek x0 ≈ π ± 1.15. Chov´an´ı algoritmu na cel´em intervalu (0, 2π) je zn´azornˇeno na obr´azku 2. Na ose X je poˇc´ateˇcn´ı podm´ınka a na ose Y je hodnota, ke kter´e algoritmus konvergoval. Pokud algoritmus divergoval3 , tak je na Y osu vynesena hodnota xk v okamˇziku n´asiln´eho ukonˇcen´ı. Je vidˇet, ˇze pro poˇca´teˇcn´ı podm´ınky bl´ızk´e 0 ˇci 2π algoritmus dokonvergoval naopak do maxima cos(x) a pro poˇca´teˇcn´ı x0 bl´ızk´e π2 ˇci 3π algoritmus odsk´akal“ kamsi do nezn´ama. 2 ” Cel´ y k´od scriptu v matlabu je v pˇr´ıloze B.
150 100
xk−>
50 0 −50 −100 −150
1
2
3
4
5
6
x0 Obr´azek 2: Chov´an´ı algoritmu Newtonovy minimalizace funkce cos(x) na intervalu poˇca´teˇcn´ıch podm´ınek (0, 2π). 2 3
. . . rozmˇeru 1x1. Pˇresnˇeji ˇreˇceno, pokud ani po 1000 iterac´ıch x1000 6= π.
3
Pˇ r´ıklad 3 Kriteri´aln´ı funkci f (x) pˇrep´ıˇsi do tvaru f (x) = xT Ax + fT x + 15, kde A =
1/2 0 0 γ/2
(8)
10 . −20
af=
Analytick´ e minimum (a) Doplnˇen´ım kvadratick´e formy (8) na ˇctverec dost´av´am analytick´e minimum 1 x = − A−1 f = 2 ∗
−10 20/γ
.
(9)
Newtonova metoda (b) Program ˇreˇs´ıc´ı tento pˇr´ıklad je v pˇr´ıloze B. Program nal´ez´a minimum v jedin´em kroku, nez´avisle na hodnotˇe γ i poˇca´teˇcn´ı podm´ınce x0 . Staˇc´ı tedy jedin´a iterace typu xk+1 = xk − H−1 gk , kde xk =
x1k x2k
(10)
, gk je gradient gk = 2Axk + f =
x1k + 10 x2k γ − 20
4
a H je Hassova matice druh´ ych derivac´ı a tedy H
−1
,
−1
= (2A)
(11) =
1 0 0 γ1
. Po dosazen´ı do
vztahu (10) dost´av´ame
1 1 1 x1 = x0 − (2A)−1 (2Ax0 + f) = x0 − A−1 2Ax0 − A−1 f = x0 − x0 − A−1 f = x∗ . 2 2 2
(12)
Vztah (12) potvrzuje, ˇze uˇz x1 , tedy odhad minima po prvn´ı iteraci Newtonova algoritmu, je pˇr´ımo minimum x∗ a to bez ohledu na parametr γ i poˇca´teˇcn´ı podm´ınku x0 .
Gradientn´ı metody (c) Gradientn´ı metoda spoˇc´ıv´a v tom, ˇze budu iterovat xk+1 = xk − αgk , α > 0, α ∈ R, 4
Pro kvadratick´ y polynom nez´ avisl´ a na bodˇe xk .
4
(13)
kde gk je gradient funkce f (x) v bodˇe xk , viz (11). V´ıcerozmˇern´a minimalizace se tak rozpad´a na opakovanou jednorozmˇernou minimalizaci pˇres parametr α. Kriteri´aln´ı funkce pro jednorozmˇernou minimalizaci, oznaˇcme ji ϕk (α), je vlastnˇe ˇrezem funkce f (x) ve smˇeru gradientu gk ϕxk (α) = f (xk − αgk ) = (xk − αgk )T A(xk − αgk ) + fT (xk − αgk ) + 15.
(14)
Nejrychlejˇ s´ı sestup (i) Pro metodu maxim´aln´ıho sestupu mus´ım urˇcit α takov´e, abych minimalizoval fci ϕk (α). Proto si vztah (14) pˇrep´ıˇsi do tvaru ϕxk (α) = α2 (gTk Agk ) + α(−gTk Axk − xTk Agk − f T gk ) + (fT xk + 15).
(15)
Doplnˇen´ım (15) na ˇctverec z´ısk´av´am, ˇze αk∗ =
gTk Axk + xTk Agk + fT gk . 2gTk Agk
(16)
V´ ysledek (16) stejn´ y, jako vztahem v pˇredn´aˇskov´ ych slidech. Jen j´a jsem k nˇemu dospˇel metodou doplnˇen´ı na ˇctverec a v pˇredn´aˇsk´ach je uk´az´ana metoda pomoc´ı derivace ϕk (α) podle α a poloˇzen´ı rovno 0. Pro ovˇeˇren´ı se d´a do pˇredn´aˇskov´eho“ vztahu ” gTk gk ∗ . (17) αk = T gk Hgk Dosadit za gk a H a po 3 aˇz 4 str´ank´ach A4 popsan´ ych hr´atky s line´arn´ı algebrou zjist´ıme rovnost vztah˚ u 5 (16) a (17). Nyn´ı k zaj´ımavˇejˇs´ı ˇca´sti - k experiment´aln´ım v´ ysledk˚ um. Na obr´azku 3 nahoˇre je vidˇet cesta“ ” algoritmu, tedy vyznaˇcen´e iterace ve vrstevnicov´em grafu. Dole na obr´azku je vynesena d´elka kroku v z´avislosti na iteraci. Je vidˇet, ˇze algoritmus chod´ı zbyteˇcnˇe klikatˇe“. ” Na dalˇs´ım obr´azku 4 je vynesen poˇcet iterac´ı, kter´e musel algoritmus vykonat, aby nalezl optimum s definovanou pˇresnost´ı. Tato pˇresnost byla stejn´a pro vˇsechny n´asleduj´ıc´ı modifikace gradientn´ıho algoritmu, aby se dali v´ ysledky navz´ajem porovn´avat. Je vidˇet, ˇze poˇcet iterac´ı je silnˇe z´avisl´ y na poˇca´teˇcn´ı podm´ınce i parametru γ, neboli ˇsiˇsatosti“ f (x). ” Program ˇreˇs´ıc´ı tento pˇr´ıklad je v pˇr´ıloze D. Backtracking (ii) Modifikace gradientn´ıho algoritmu pomoc´ı backtrackingu spoˇc´ıv´a jen v jin´em zp˚ usobu nalezen´ı d´elky kroku - koeficientu α. Naprogramoval jsem hled´an´ı α pˇresnˇe tak, jak je uvedeno v pˇredn´aˇskov´ ych slidech a nahradil t´ım v´ ypoˇcet podle vztahu (16). V´ ysledky jsou na obr´azc´ıch 5 a 6. Na rozd´ıl od pˇredchoz´ıho pˇr´ıpadu trv´a algoritmus pr˚ umˇernˇe v´ıc iterac´ı, ale odstranila se z´avislost na poˇca´teˇcn´ı podm´ınce a ˇca´steˇcnˇe i na parametru γ. Program ˇreˇs´ıc´ı tento pˇr´ıklad je v pˇr´ıloze E. 5
Za pˇredpokladu symetrick´e matice H! Tedy HT = H . . .
5
gamma = 10
gamma = 30 10
10
5 x2
x2
5
0
0
−5
−5 −15
−10 x1
−5
0
−20
4
8
3
6
delka kroku
delka kroku
−20
2 1 0
0
20
40 iterace
60
−10 x1
−5
0
4 2 0
80
−15
0
10
20
30
iterace
Obr´azek 3: Gradientn´ı metoda minimalizace f (x) pomoc´ı nejrychlejˇs´ıho sestupu. Pr˚ ubˇeh iterac´ı pro ϕ = 20 (nahoˇre) a z´avislost d´elky kroku na iteraci (dole). gamma = 10
gamma = 30 300
80
pocet iteraci
pocet iteraci
100
60 40
200
100
20 0
0
100
200 fi [deg]
300
0
400
0
100
200 fi [deg]
300
400
Obr´azek 4: Poˇcet iterac´ı algoritmu gradientn´ı minimalizace f (x) pomoc´ı nejrychlejˇs´ıho sestupu v z´avislosti na ϕ. PARTAN (iii) Pˇri u ´pravˇe gradientn´ıho algoritmu metodou PARTAN jsem vych´azel z metody nejrychlejˇs´ıho sestupu. Stˇr´ıd´a se krok klasick´eho algoritmu s krokem zrychlen´ ym. Protoˇze ke kroku zrychlen´emu je potˇreba m´ıt alespoˇ n dva pˇredchoz´ı kroky, na zaˇca´tku algoritmu se v´ yjimeˇcnˇe provedou dva 6
gamma = 10
gamma = 30 10
10
5 x2
x2
5
0
0
−5
−5 −20
−15
−10 x1
−5
0
−20
4
2
0
−10 x1
−5
0
3 delka kroku
delka kroku
6
−15
0
20
40
2
1
0
60
0
50 iterace
iterace
100
Obr´azek 5: Gradientn´ı metoda minimalizace f (x) pomoc´ı backtrackingu. Pr˚ ubˇeh iterac´ı pro ϕ = 20 (nahoˇre) a z´avislost d´elky kroku na iteraci (dole). gamma = 10
gamma = 30 100 80
60
pocet iteraci
pocet iteraci
80
40 20 0
60 40 20
0
100
200 fi [deg]
300
0
400
0
100
200 fi [deg]
300
400
Obr´azek 6: Poˇcet iterac´ı algoritmu gradientn´ı minimalizace f (x) pomoc´ı backtrackingu v z´avislosti na ϕ. klasick´e kroky. V´ ysledky jsou na obr´azc´ıch 7 a 8. V tomto pˇr´ıpadˇe pˇrin´aˇs´ı metoda obrovsk´e zlepˇsen´ı, protoˇze souˇcet smˇer˚ u dvou pˇredchoz´ıch gradient˚ u d´av´a smˇer ukazuj´ıc´ı pˇr´ımo do minima. Protoˇze i pˇri ∗ zrychlen´em kroku se hled´a αk metodou nejvˇetˇs´ıho sestupu, viz. (16), dos´ahneme 3. krokem vˇzdy 7
optima. gamma = 10
gamma = 30 10
10
5 x2
x2
5 0
0 −5
−5 −20
−15
−10 x1
−5
0
−20
−15
−10 x1
−5
0
Obr´azek 7: Pr˚ ubˇeh iterac´ı gradientn´ıho algoritmu minimalizace f (x) pomoc´ı metody PARTAN pro ϕ = 20. gamma = 30
4
4
3.5
3.5
pocet iteraci
pocet iteraci
gamma = 10
3 2.5 2
0
100
200 fi [deg]
300
3 2.5 2
400
0
100
200 fi [deg]
300
400
Obr´azek 8: Poˇcet iterac´ı gradientn´ıho algoritmu minimalizace f (x) pomoc´ı metody PARTAN v z´avislosti na ϕ. Program ˇreˇs´ıc´ı tento pˇr´ıklad je v pˇr´ıloze F.
Minimalizace pomoc´ı matlabu (iiii) Vyˇreˇsil jsem minimalizaci f (x) pomoc´ı matlabsk´e funkce fminsearch. Tato funkce pouˇz´ıv´a algoritmus flexibiln´ıho simplexu oznaˇcovan´eho t´eˇz Nelder-Mead ˇci v´ ystiˇznˇe Am´eba. V´ yhoda tohoto algoritmu je, ˇze nepotˇrebuje zn´at derivace funkce. Proto je jeˇstˇe obecnˇejˇs´ı, neˇz gradientn´ı metody. V´ ysledky jsou vidˇet na obr´azc´ıch 9 a 10. Algoritmus mˇe pˇrekvapil relativnˇe mal´ ym poˇctem iterac´ı, kter´e potˇrebuje k nalezen´ı optima. Jedn´a se o srovnateln´e ˇc´ıslo, jako u gradientn´ıch metod. Algoritmus ale nen´ı pˇr´ıliˇs bezpeˇcn´ y“. Pro urˇcit´e poˇca´teˇcn´ı podm´ınky se stalo, ˇze poˇcet ” iterac´ı potˇrebn´ ych k nalezen´ı minima doslova explodoval“. V jednom pˇr´ıpadˇe bylo zapotˇreb´ı ” 12000 iterac´ı. 8
gamma = 10
gamma = 30 10
10
5 x2
x2
5 0
0 −5
−5 −20
−15
−10 x1
−5
0
−20
−15
−10 x1
−5
0
Obr´azek 9: Kroky minimalizaˇcn´ıho algoritmu Am´eba na funkci f (x). gamma = 10
gamma = 30 500 pocet iteraci
pocet iteraci
400 300 200 100 0
400 300 200 100
0
100
200 fi [deg]
300
0
400
0
100
200 fi [deg]
300
Obr´azek 10: Poˇcet iterac´ı algoritmu Am´eba v z´avislosti na ϕ. Program ˇreˇs´ıc´ı tento pˇr´ıklad je v pˇr´ıloze G.
9
400
Pˇ r´ılohy Pˇ r´ıloha A Program v PHP pro v´ ypoˇcet
√ 13 na N desetinn´ ych m´ıst.
\n"; } //vysledne vypsani cisla sqrt(13) echo $xk; ?>
Pˇ r´ıloha B Script v matlabu ˇreˇs´ıc´ı pˇr´ıklad 2a. clear; clc; close all; x0 = pi-2; a = [0:0.01:2*pi]; %% vykresleni grafu cos(x) a aproximace pomoci g(x) v bode x0 10
plot(a,cos(a)); hold on; grid on; g = cos(x0) - sin(x0)*(a-x0) - cos(x0)*(a-x0).^2./2; plot(a,g,’red’); axis([0 2*pi -1.5 1.5]); xk = x0 - sin(x0)/cos(x0); %jeden krok newtonova alg. hold on; plot([pi pi], [0 -1],’b --’); plot([x0 x0], [0 cos(x0)],’k --’); plot([xk xk], [0 cos(x0) - sin(x0)*(xk-x0) - cos(x0)*(xk-x0)^2/2;],’r --’); legend(’cos(x)’,’g(x)’,’x*’,’xk’,’xk+1’) xlabel(’x’); %% hledani minima cos(x) pro ruzne poc. podminky x0 x0s = [0.01:0.01:2*pi-0.01]; xopts = x0s; for k = 1:size(x0s,2) n = 0; while (abs(xopts(k) - pi) > 0.00001) && n < 1000 xopts(k) = xopts(k) - sin(xopts(k))/cos(xopts(k)); n = n + 1; end end figure(2); plot(x0s,xopts); grid on; xlabel(’x0’); ylabel(’xk->’);
Pˇ r´ıloha C Script v matlabu ˇreˇs´ıc´ı pˇr´ıklad 3b. clear; clc; close all; s = 30; %gamma A = [1/2 0; 0 s/2]; f = [10; -20]; xopt = [-10; 20/s]; r = 10; fi = pi*20/180; x0 = [r*cos(fi)+xopt(1); r*sin(fi)+xopt(2)]; iteraci = 0; xk = x0; H = [1 0; 0 s]; 11
while abs(xk - xopt) > eps*ones(2,1) gk = [xk(1)+10; s*xk(2)-20]; xk = xk - inv(H)*gk; iteraci = iteraci + 1; end disp([’Optimum after ’ num2str(iteraci) ’ iterrations is x=[’ num2str(xk(1)) ’, ’ num2str(xk(2)) ’].’]);
Pˇ r´ıloha D Script v matlabu ˇreˇs´ıc´ı pˇr´ıklad 3ci . clear; clc; close all; f = [10;-20]; r = 10; fi = pi*10/180; presnost = 0.00001; %% vykresleni vrstevnicovych grafu % gamma = 10 s = 10; xopt = [-10; 20/s]; x1 = [(xopt(1)-10):0.05:(xopt(1)+10)]; x2 = [(xopt(2)-10):0.05:(xopt(2)+10)]; [X1,X2] = meshgrid(x1, x2); J = zeros(size(X1)); for a = 1 : size(X1,1) for b = 1 : size(X1,2) x1i = X1(a,b); x2i = X2(a,b); J(a,b) = (x1i^2+s*x2i^2)/2+10*x1i-20*x2i+15; end; end; subplot(2,2,1); contour(X1,X2,J); xlabel (’x1’); ylabel (’x2’); title(’gamma = 10’); % gamma = 30 s = 30; xopt = [-10; 20/s]; 12
x1 = [(xopt(1)-10):0.05:(xopt(1)+10)]; x2 = [(xopt(2)-10):0.05:(xopt(2)+10)]; [X1,X2] = meshgrid(x1, x2); J = zeros(size(X1)); for a = 1 : size(X1,1) for b = 1 : size(X1,2) x1i = X1(a,b); x2i = X2(a,b); J(a,b) = (x1i^2+s*x2i^2)/2+10*x1i-20*x2i+15; end; end; subplot(2,2,2); contour(X1,X2,J); xlabel (’x1’); ylabel (’x2’); title(’gamma = 30’); %% nalezeni optima pomoci gradinetni metody nejvetsiho sestupu % gamma = 10 s = 10; xopt = [-10; 20/s]; A = [1/2 0; 0 s/2]; x0 = [r*cos(fi)+xopt(1); r*sin(fi)+xopt(2)]; iteraci = 0; cesta = x0; delky_kroku = 0; xk = x0; tolerance = presnost*2; while tolerance > presnost gk = [xk(1) + 10; s*xk(2) - 20]; %urcim smer = -gradient a = (gk’*A*xk + xk’*A*gk + f’*gk) / (2*gk’*A*gk); %urcim optimalni delku kroku xk = xk - a*gk; iteraci = iteraci+1; cesta(:,iteraci+1) = xk; delky_kroku(iteraci+1) = sqrt(sum((cesta(:,iteraci)-xk).^2)); tolerance = abs((xopt’*A*xopt - f’*xopt - 15)-(xk’*A*xk - f’*xk - 15)); end disp([’Optimum with tolerance ’ num2str(tolerance) ’ founded after ’ num2str(iteraci) ’ figure(1); subplot(2,2,1); hold on; plot(cesta(1,:), cesta(2,:),’k+-’); 13
subplot(2,2,3); plot([1:iteraci+1], delky_kroku); xlabel(’iterace’); ylabel(’delka kroku’); grid on; %vykresleni poctu iteraci v zavislosti na fi figure(2); subplot(1,2,1); fis = 1:360; poc_i = [fis; zeros(size(fis))]; for fi = fis x0 = [r*cos(pi*fi/180)+xopt(1); r*sin(pi*fi/180)+xopt(2)]; iteraci = 0; xk = x0; tolerance = presnost*2; while tolerance > presnost gk = [xk(1) + 10; s*xk(2) - 20]; %urcim smer = -gradient a = (gk’*A*xk + xk’*A*gk + f’*gk) / (2*gk’*A*gk); %urcim optimalni delku kroku xk = xk - a*gk; iteraci = iteraci+1; tolerance = abs((xopt’*A*xopt - f’*xopt - 15)-(xk’*A*xk - f’*xk - 15)); end poc_i(2,fi) = iteraci; end plot(poc_i(1,:), poc_i(2,:)); xlabel(’fi [deg]’); ylabel(’pocet iteraci’); grid on; title(’gamma = 10’); % gamma = 30 s = 30; fi = pi*10/180; A = [1/2 0; 0 s/2]; x0 = [r*cos(fi)+xopt(1); r*sin(fi)+xopt(2)]; xopt = [-10; 20/s]; iteraci = 0; cesta = x0; delky_kroku = 0; xk = x0; tolerance = presnost*2; while tolerance > presnost gk = [xk(1) + 10; s*xk(2) - 20]; %urcim smer = -gradient a = (gk’*A*xk + xk’*A*gk + f’*gk) / (2*gk’*A*gk); %urcim optimalni delku kroku xk = xk - a*gk; 14
iteraci = iteraci+1; cesta(:,iteraci+1) = xk; delky_kroku(iteraci+1) = sqrt(sum((cesta(:,iteraci)-xk).^2)); tolerance = abs((xopt’*A*xopt - f’*xopt - 15)-(xk’*A*xk - f’*xk - 15)); end disp([’Optimum with tolerance ’ num2str(tolerance) ’ founded after ’ num2str(iteraci) ’ figure(1); subplot(2,2,2); hold on; plot(cesta(1,:), cesta(2,:),’k+-’); subplot(2,2,4); plot([1:iteraci+1], delky_kroku); xlabel(’iterace’); ylabel(’delka kroku’); grid on; %vykresleni poctu iteraci v zavislosti na fi figure(2); subplot(1,2,2); fis = 1:360; poc_i = [fis; zeros(size(fis))]; for fi = fis x0 = [r*cos(pi*fi/180)+xopt(1); r*sin(pi*fi/180)+xopt(2)]; iteraci = 0; xk = x0; tolerance = presnost*2; while tolerance > presnost gk = [xk(1) + 10; s*xk(2) - 20]; %urcim smer = -gradient a = (gk’*A*xk + xk’*A*gk + f’*gk) / (2*gk’*A*gk); %urcim optimalni delku kroku xk = xk - a*gk; iteraci = iteraci+1; tolerance = abs((xopt’*A*xopt - f’*xopt - 15)-(xk’*A*xk - f’*xk - 15)); end poc_i(2,fi) = iteraci; end plot(poc_i(1,:), poc_i(2,:)); xlabel(’fi [deg]’); ylabel(’pocet iteraci’); grid on; title(’gamma = 30’);
Pˇ r´ıloha E Script v matlabu ˇreˇs´ıc´ı pˇr´ıklad 3cii .
15
clear; clc; close all; f = [10;-20]; r = 10; fi = pi*10/180; presnost = 0.00001; %% vykresleni vrstevnicovych grafu % gamma = 10 s = 10; xopt = [-10; 20/s]; x1 = [(xopt(1)-10):0.05:(xopt(1)+10)]; x2 = [(xopt(2)-10):0.05:(xopt(2)+10)]; [X1,X2] = meshgrid(x1, x2); J = zeros(size(X1)); for a = 1 : size(X1,1) for b = 1 : size(X1,2) x1i = X1(a,b); x2i = X2(a,b); J(a,b) = (x1i^2+s*x2i^2)/2+10*x1i-20*x2i+15; end; end; subplot(2,2,1); contour(X1,X2,J); xlabel (’x1’); ylabel (’x2’); title(’gamma = 10’); % gamma = 30 s = 30; xopt = [-10; 20/s]; x1 = [(xopt(1)-10):0.05:(xopt(1)+10)]; x2 = [(xopt(2)-10):0.05:(xopt(2)+10)]; [X1,X2] = meshgrid(x1, x2); J = zeros(size(X1)); for a = 1 : size(X1,1) for b = 1 : size(X1,2) x1i = X1(a,b); x2i = X2(a,b); J(a,b) = (x1i^2+s*x2i^2)/2+10*x1i-20*x2i+15; end; end;
16
subplot(2,2,2); contour(X1,X2,J); xlabel (’x1’); ylabel (’x2’); title(’gamma = 30’); %% nalezeni optima pomoci gradinetni metody + backtrackingu % gamma = 10 s = 10; xopt = [-10; 20/s]; A = [1/2 0; 0 s/2]; x0 = [r*cos(fi)+xopt(1); r*sin(fi)+xopt(2)];
iteraci = 0; cesta = x0; delky_kroku = 0; xk = x0; tolerance = presnost*2; while tolerance > presnost gk = [xk(1) + 10; s*xk(2) - 20]; %urcim smer = -gradient %hledani a pomoci backtracking e = 0.2; b = 0.5; a = 1; while ([xk-a*gk]’*A*[xk-a*gk] + f’*[xk-a*gk] + 15) > ((xk’*A*xk + f’*xk + 15) - e*a* a = b*a; end xk = xk - a*gk; iteraci = iteraci+1; cesta(:,iteraci+1) = xk; delky_kroku(iteraci+1) = sqrt(sum((cesta(:,iteraci)-xk).^2)); tolerance = abs((xopt’*A*xopt - f’*xopt - 15)-(xk’*A*xk - f’*xk - 15)); end disp([’Optimum with tolerance ’ num2str(tolerance) ’ founded after ’ num2str(iteraci) ’ figure(1); subplot(2,2,1); hold on; plot(cesta(1,:), cesta(2,:),’k+-’); subplot(2,2,3); plot([1:iteraci+1], delky_kroku); xlabel(’iterace’); ylabel(’delka kroku’); grid on; %vykresleni poctu iteraci v zavislosti na fi figure(2); 17
subplot(1,2,1); fis = 1:360; poc_i = [fis; zeros(size(fis))]; for fi = fis x0 = [r*cos(pi*fi/180)+xopt(1); r*sin(pi*fi/180)+xopt(2)]; iteraci = 0; xk = x0; tolerance = presnost*2; while tolerance > presnost gk = [xk(1) + 10; s*xk(2) - 20]; %urcim smer = -gradient %hledani a pomoci backtracking e = 0.2; b = 0.5; a = 1; while ([xk-a*gk]’*A*[xk-a*gk] + f’*[xk-a*gk] + 15) > ((xk’*A*xk + f’*xk + 15) a = b*a; end xk = xk - a*gk; iteraci = iteraci+1; tolerance = abs((xopt’*A*xopt - f’*xopt - 15)-(xk’*A*xk - f’*xk - 15)); end poc_i(2,fi) = iteraci; end plot(poc_i(1,:), poc_i(2,:)); xlabel(’fi [deg]’); ylabel(’pocet iteraci’); grid on; title(’gamma = 10’); % gamma = 30 s = 30; fi = pi*10/180; A = [1/2 0; 0 s/2]; x0 = [r*cos(fi)+xopt(1); r*sin(fi)+xopt(2)]; xopt = [-10; 20/s]; iteraci = 0; cesta = x0; delky_kroku = 0; xk = x0; tolerance = presnost*2;
while tolerance > presnost gk = [xk(1) + 10; s*xk(2) - 20]; %urcim smer = -gradient %hledani a pomoci backtracking e = 0.2; b = 0.5; a = 1; while ([xk-a*gk]’*A*[xk-a*gk] + f’*[xk-a*gk] + 15) > ((xk’*A*xk + f’*xk + 15) - e*a* 18
a = b*a; end xk = xk - a*gk; iteraci = iteraci+1; cesta(:,iteraci+1) = xk; delky_kroku(iteraci+1) = sqrt(sum((cesta(:,iteraci)-xk).^2)); tolerance = abs((xopt’*A*xopt - f’*xopt - 15)-(xk’*A*xk - f’*xk - 15)); end disp([’Optimum with tolerance ’ num2str(tolerance) ’ founded after ’ num2str(iteraci) ’ figure(1); subplot(2,2,2); hold on; plot(cesta(1,:), cesta(2,:),’k+-’); subplot(2,2,4); plot([1:iteraci+1], delky_kroku); xlabel(’iterace’); ylabel(’delka kroku’); grid on; %vykresleni poctu iteraci v zavislosti na fi figure(2); subplot(1,2,2); fis = 1:360; poc_i = [fis; zeros(size(fis))]; for fi = fis x0 = [r*cos(pi*fi/180)+xopt(1); r*sin(pi*fi/180)+xopt(2)]; iteraci = 0; xk = x0; tolerance = presnost*2; while tolerance > presnost gk = [xk(1) + 10; s*xk(2) - 20]; %urcim smer = -gradient %hledani a pomoci backtracking e = 0.2; b = 0.5; a = 1; while ([xk-a*gk]’*A*[xk-a*gk] + f’*[xk-a*gk] + 15) > ((xk’*A*xk + f’*xk + 15) a = b*a; end xk = xk - a*gk; iteraci = iteraci+1; tolerance = abs((xopt’*A*xopt - f’*xopt - 15)-(xk’*A*xk - f’*xk - 15)); end poc_i(2,fi) = iteraci; end plot(poc_i(1,:), poc_i(2,:)); 19
xlabel(’fi [deg]’); ylabel(’pocet iteraci’); grid on; title(’gamma = 30’);
Pˇ r´ıloha F Script v matlabu ˇreˇs´ıc´ı pˇr´ıklad 3ciii . clear; clc; close all; f = [10;-20]; r = 10; fi = pi*10/180; presnost = 0.00001; %% vykresleni vrstevnicovych grafu % gamma = 10 s = 10; xopt = [-10; 20/s]; x1 = [(xopt(1)-10):0.05:(xopt(1)+10)]; x2 = [(xopt(2)-10):0.05:(xopt(2)+10)]; [X1,X2] = meshgrid(x1, x2); J = zeros(size(X1)); for a = 1 : size(X1,1) for b = 1 : size(X1,2) x1i = X1(a,b); x2i = X2(a,b); J(a,b) = (x1i^2+s*x2i^2)/2+10*x1i-20*x2i+15; end; end; subplot(1,2,1); contour(X1,X2,J); xlabel (’x1’); ylabel (’x2’); title(’gamma = 10’); % gamma = 30 s = 30; xopt = [-10; 20/s]; x1 = [(xopt(1)-10):0.05:(xopt(1)+10)]; x2 = [(xopt(2)-10):0.05:(xopt(2)+10)]; [X1,X2] = meshgrid(x1, x2); 20
J = zeros(size(X1)); for a = 1 : size(X1,1) for b = 1 : size(X1,2) x1i = X1(a,b); x2i = X2(a,b); J(a,b) = (x1i^2+s*x2i^2)/2+10*x1i-20*x2i+15; end; end; subplot(1,2,2); contour(X1,X2,J); xlabel (’x1’); ylabel (’x2’); title(’gamma = 30’); %% nalezeni optima pomoci gradinetni metody + backtrackingu % gamma = 10 s = 10; xopt = [-10; 20/s]; A = [1/2 0; 0 s/2]; x0 = [r*cos(fi)+xopt(1); r*sin(fi)+xopt(2)]; iteraci = 0; cesta = x0; delky_kroku = 0; xk = x0; tolerance = presnost*2; %prvni iterace xk_old = xk; gk = [xk(1) + 10; s*xk(2) - 20]; %urcim smer = -gradient a = (gk’*A*xk + xk’*A*gk + f’*gk) / (2*gk’*A*gk); %urcim optimalni delku kroku xk = xk - a*gk; iteraci = iteraci+1; cesta(:,iteraci+1) = xk; while tolerance > presnost gk = [xk(1) + 10; s*xk(2) - 20]; %urcim smer = -gradient %nejmensi sestup a = (gk’*A*xk + xk’*A*gk + f’*gk) / (2*gk’*A*gk); %urcim optimalni delku kroku yk = xk - a*gk;%pomocny bod iteraci = iteraci+1; cesta(:,iteraci+1) = yk; gk = xk_old-yk; xk_old = xk; 21
%nejmensi sestup - druhe hledani a = (gk’*A*yk + yk’*A*gk + f’*gk) / (2*gk’*A*gk); %urcim optimalni delku kroku xk = yk - a*gk; %zrychleny krok iteraci = iteraci+1; cesta(:,iteraci+1) = xk; tolerance = abs((xopt’*A*xopt - f’*xopt - 15)-(xk’*A*xk - f’*xk - 15)); end disp([’Optimum with tolerance ’ num2str(tolerance) ’ founded after ’ num2str(iteraci) ’ figure(1); subplot(1,2,1); hold on; plot(cesta(1,:), cesta(2,:),’k+-’);
%vykresleni poctu iteraci v zavislosti na fi figure(2); subplot(1,2,1); fis = 1:360; poc_i = [fis; zeros(size(fis))]; for fi = fis x0 = [r*cos(pi*fi/180)+xopt(1); r*sin(pi*fi/180)+xopt(2)]; iteraci = 0; xk = x0; tolerance = presnost*2; %prvni iterace xk_old = xk; gk = [xk(1) + 10; s*xk(2) - 20]; %urcim smer = -gradient a = (gk’*A*xk + xk’*A*gk + f’*gk) / (2*gk’*A*gk); %urcim optimalni delku kroku xk = xk - a*gk; iteraci = iteraci+1; cesta(:,iteraci+1) = xk; while tolerance > presnost gk = [xk(1) + 10; s*xk(2) - 20]; %urcim smer = -gradient %nejmensi sestup a = (gk’*A*xk + xk’*A*gk + f’*gk) / (2*gk’*A*gk); %urcim optimalni delku kroku yk = xk - a*gk;%pomocny bod iteraci = iteraci+1; cesta(:,iteraci+1) = yk; gk = xk_old-yk; 22
xk_old = xk; %nejmensi sestup - druhe hledani a = (gk’*A*yk + yk’*A*gk + f’*gk) / (2*gk’*A*gk); %urcim optimalni delku kroku xk = yk - a*gk; %zrychleny krok iteraci = iteraci+1; cesta(:,iteraci+1) = xk; tolerance = abs((xopt’*A*xopt - f’*xopt - 15)-(xk’*A*xk - f’*xk - 15)); end
poc_i(2,fi) = iteraci; end plot(poc_i(1,:), poc_i(2,:)); xlabel(’fi [deg]’); ylabel(’pocet iteraci’); grid on; title(’gamma = 10’); % gamma = 30 s = 30; fi = pi*10/180; xopt = [-10; 20/s]; A = [1/2 0; 0 s/2]; x0 = [r*cos(fi)+xopt(1); r*sin(fi)+xopt(2)]; iteraci = 0; cesta = x0; delky_kroku = 0; xk = x0; tolerance = presnost*2; %prvni iterace xk_old = xk; gk = [xk(1) + 10; s*xk(2) - 20]; %urcim smer = -gradient a = (gk’*A*xk + xk’*A*gk + f’*gk) / (2*gk’*A*gk); %urcim optimalni delku kroku xk = xk - a*gk; iteraci = iteraci+1; cesta(:,iteraci+1) = xk; while tolerance > presnost gk = [xk(1) + 10; s*xk(2) - 20]; %urcim smer = -gradient %nejmensi sestup a = (gk’*A*xk + xk’*A*gk + f’*gk) / (2*gk’*A*gk); %urcim optimalni delku kroku yk = xk - a*gk;%pomocny bod iteraci = iteraci+1; cesta(:,iteraci+1) = yk; 23
gk = xk_old-yk; xk_old = xk; %nejmensi sestup - druhe hledani a = (gk’*A*yk + yk’*A*gk + f’*gk) / (2*gk’*A*gk); %urcim optimalni delku kroku xk = yk - a*gk; %zrychleny krok iteraci = iteraci+1; cesta(:,iteraci+1) = xk; tolerance = abs((xopt’*A*xopt - f’*xopt - 15)-(xk’*A*xk - f’*xk - 15)); end disp([’Optimum with tolerance ’ num2str(tolerance) ’ founded after ’ num2str(iteraci) ’ figure(1); subplot(1,2,2); hold on; plot(cesta(1,:), cesta(2,:),’k+-’);
%vykresleni poctu iteraci v zavislosti na fi figure(2); subplot(1,2,2); fis = 1:360; poc_i = [fis; zeros(size(fis))]; for fi = fis x0 = [r*cos(pi*fi/180)+xopt(1); r*sin(pi*fi/180)+xopt(2)]; iteraci = 0; xk = x0; tolerance = presnost*2; %prvni iterace xk_old = xk; gk = [xk(1) + 10; s*xk(2) - 20]; %urcim smer = -gradient a = (gk’*A*xk + xk’*A*gk + f’*gk) / (2*gk’*A*gk); %urcim optimalni delku kroku xk = xk - a*gk; iteraci = iteraci+1; cesta(:,iteraci+1) = xk; while tolerance > presnost gk = [xk(1) + 10; s*xk(2) - 20]; %urcim smer = -gradient %nejmensi sestup a = (gk’*A*xk + xk’*A*gk + f’*gk) / (2*gk’*A*gk); %urcim optimalni delku kroku yk = xk - a*gk;%pomocny bod iteraci = iteraci+1; 24
cesta(:,iteraci+1) = yk; gk = xk_old-yk; xk_old = xk; %nejmensi sestup - druhe hledani a = (gk’*A*yk + yk’*A*gk + f’*gk) / (2*gk’*A*gk); %urcim optimalni delku kroku xk = yk - a*gk; %zrychleny krok iteraci = iteraci+1; cesta(:,iteraci+1) = xk; tolerance = abs((xopt’*A*xopt - f’*xopt - 15)-(xk’*A*xk - f’*xk - 15)); end
poc_i(2,fi) = iteraci; end plot(poc_i(1,:), poc_i(2,:)); xlabel(’fi [deg]’); ylabel(’pocet iteraci’); grid on; title(’gamma = 30’);
Pˇ r´ıloha G Script v matlabu ˇreˇs´ıc´ı pˇr´ıklad 3ciiii . clear; clc; close all; f = [10;-20]; r = 10; fi = pi*10/180; presnost = 0.00001; %% vykresleni vrstevnicovych grafu % gamma = 10 s = 10; xopt = [-10; 20/s]; x1 = [(xopt(1)-10):0.05:(xopt(1)+10)]; x2 = [(xopt(2)-10):0.05:(xopt(2)+10)]; [X1,X2] = meshgrid(x1, x2); J = zeros(size(X1)); for a = 1 : for b = x1i x2i
size(X1,1) 1 : size(X1,2) = X1(a,b); = X2(a,b); 25
J(a,b) = (x1i^2+s*x2i^2)/2+10*x1i-20*x2i+15; end; end; subplot(1,2,1); contour(X1,X2,J); xlabel (’x1’); ylabel (’x2’); title(’gamma = 10’); % gamma = 30 s = 30; xopt = [-10; 20/s]; x1 = [(xopt(1)-10):0.05:(xopt(1)+10)]; x2 = [(xopt(2)-10):0.05:(xopt(2)+10)]; [X1,X2] = meshgrid(x1, x2); J = zeros(size(X1)); for a = 1 : size(X1,1) for b = 1 : size(X1,2) x1i = X1(a,b); x2i = X2(a,b); J(a,b) = (x1i^2+s*x2i^2)/2+10*x1i-20*x2i+15; end; end; subplot(1,2,2); contour(X1,X2,J); xlabel (’x1’); ylabel (’x2’); title(’gamma = 30’); %% nalezeni optima pomoci ameby % gamma = 10 s = 10; xopt = [-10; 20/s]; A = [1/2 0; 0 s/2]; x0 = [r*cos(fi)+xopt(1); r*sin(fi)+xopt(2)];
global cesta cesta = [x0]; optiony = optimset(’TolFun’,presnost, ’MaxFunEvals’, 100000, ’MaxIter’, 100000); %vlastni algoritmus hledani xk = fminsearch(@(x) pr3fun(x,s),x0,optiony); 26
tolerance = abs((xopt’*A*xopt - f’*xopt - 15)-(xk’*A*xk - f’*xk - 15)); disp([’Optimum with tolerance ’ num2str(tolerance) ’ founded after ’ num2str(size(cesta, figure(1); subplot(1,2,1); hold on; plot(cesta(1,:), cesta(2,:),’k+-’); %vykresleni poctu iteraci v zavislosti na fi figure(2); subplot(1,2,1); fis = 1:360; poc_i = [fis; zeros(size(fis))]; for fi = fis x0 = [r*cos(pi*fi/180)+xopt(1); r*sin(pi*fi/180)+xopt(2)]; global cesta cesta = [x0]; %vlastni algoritmus hledani xk = fminsearch(@(x) pr3fun(x,s),x0,optiony); poc_i(2,fi) = size(cesta,2); end plot(poc_i(1,:), poc_i(2,:)); xlabel(’fi [deg]’); ylabel(’pocet iteraci’); grid on; title(’gamma = 10’); % gamma = 30 s = 30; fi = pi*10/180; xopt = [-10; 20/s]; A = [1/2 0; 0 s/2]; x0 = [r*cos(fi)+xopt(1); r*sin(fi)+xopt(2)];
global cesta cesta = [x0]; %vlastni algoritmus hledani xk = fminsearch(@(x) pr3fun(x,s),x0,optiony);
tolerance = abs((xopt’*A*xopt - f’*xopt - 15)-(xk’*A*xk - f’*xk - 15)); disp([’Optimum with tolerance ’ num2str(tolerance) ’ founded after ’ num2str(size(cesta, figure(1); subplot(1,2,2); 27
hold on; plot(cesta(1,:), cesta(2,:),’k+-’); %vykresleni poctu iteraci v zavislosti na fi figure(2); subplot(1,2,2); fis = 1:360; poc_i = [fis; zeros(size(fis))]; for fi = fis x0 = [r*cos(pi*fi/180)+xopt(1); r*sin(pi*fi/180)+xopt(2)]; global cesta cesta = [x0]; %vlastni algoritmus hledani xk = fminsearch(@(x) pr3fun(x,s),x0,optiony); poc_i(2,fi) = size(cesta,2); end plot(poc_i(1,:), poc_i(2,:)); xlabel(’fi [deg]’); ylabel(’pocet iteraci’); grid on; title(’gamma = 30’);
28
