Post on 24-Dec-2019
transcript
5
Lucrare de laborator Minimizarea funcţiilor de o singură variabilă
în absenţa restricţiilor
1. Scopul lucrării
Principalul obiectiv al acestei lucrări este de a prezenta modul de
lansare a unei probleme de minimizare a unei funcţii obiectiv de o singură
variabilă în absenţa restricţiilor precum şi soluţionarea unei astfel de
probleme utilizând tehnica de calcul numerică. Prezentarea este propusă în
cadrul pachetului de programe Matlab (versiunea 6.5). De altfel, un al
doilea obiectiv al acestei lucrări este familiarizarea utilizatorului cu
programul amintit care reprezintă probabil cel mai elaborat program pentru
calcule ştiinţifice.
2. Prezentarea lucrării
Problema minimizării unei funcţii de o singură variabilă în absenţa
restricţiilor se formulează simplu în forma: 1
min argmin ( ) , [ , ]x f x x a b R= ∈ ⊂ .
Problema, simplă la prima vedere, poate fi foarte complicată dacă
funcţia obiectiv are mai multe limite locale. Într-un astfel de caz majoritatea
algoritmilor de căutare a minimului fie se blochează, fie converg către o
soluţie de minim local ceea ce nu oferă o soluţie admisibilă. În acest sens, în
cele ce urmează vom considera că intervalul [ ]ba, constituie un interval de
incertitudine ce reprezintă un interval în care se găseşte cu certitudine un
minim unic al funcţiei obiectiv. Sigur, o astfel de condiţionare este
restrictivă, dar în majoritatea aplicaţiilor putem estima un asemenea interval.
1
6
Principalele subrutine conţinute în Optimization Toolbox din cadrul
pachetului de programe Matlab ce permit soluţionarea problemei propuse
sunt:
• fminbnd – pentru probleme de optimizare scalară;
• fminsearch – pentru probleme de optimizare neliniară în absenţa
restricţiilor.
În continuare, vom prezenta pe scurt sintaxa cu care aceste subrutine sunt
apelate.
Apelarea primei subrutine, în funcţie de necesităţi, poate fi făcută în
una din următoarele variante:
1) ( )21,,min xxfunbndfx =
2) ( )optionsxxfunbndfx ,,,min 21=
3) ( ),...,,,,,min 2121 PPoptionsxxfunbndfx =
Prima formă de apelare întoarce c rezultat minimul funcţiei fun din
intervalul [ ]21 , xx . Cea de a doua formă de apelare permite evaluarea
aceluiaşi minim determinat în condiţiile impuse prin structura options.
Ultima formă de apelare se impune în cazul în care funcţia conţine mai
mulţi parametrii.
Opţiunile ce pot fi impuse pot fi:
• Display – fixează nivelul prezentării. Acesta poate fi: off în cazul în care
nu dorim informaţii suplimentare, iter în cazul în care urmărim evoluţia
procesului iterativ, final caz în care este prezentată numai evoluţia finală
• TolX prin care se impune toleranţa finală admisă;
• MaxFunEvals – se impune numărul maxim de funcţii ce pot fi evaluate;
• MaxIter – număr maxim de iteraţii.
7
Impunerea structurii Options se realizează astfel:
options=optimset(‘param1’,value1,’param2’,value2,…).
Subrutina fminsearch va fi apelată în una din formele:
1) ( )0,min xfunsearchf ;
2) ( )optionsxfunsearchf ,,min 0 ;
3) ( ),...2,1,,,min 0 PPoptionsxfunsearchf .
Variabila de intrare 0x reprezintă punctul de iniţializare a căutării.
Celelalte variabile de intrare au semnificaţia anterior prezentată.
În continuare vom prezenta un exemplu simplu privitor la problema
propusă. Considerăm funcţia obiectiv:
( ) ( )xxxf −⋅−= 1 .
Se caută punctul de minim situat în intervalul [ ]10,0 . Pentru început,
realizăm o m-funcţie asociată funcţiei impuse pe care o definim f1.
function y=f1(x)
y=-x.*(1-x);
În continuare impunem structura options:
>>options=optimset('TolX',1e-12,'Display','iter','TolFun',1e-8);
Lansăm procedura de căutare:
>> x=fminbnd('f1',0,10)
Rezultatele afişate sunt prezentate în cele ce urmează: Func-count x f(x) Procedure
1 3.81966 10.7701 initial
2 6.18034 32.0163 golden
3 2.36068 3.21213 golden
4 0.5 -0.25 parabolic
8
5 0.5 -0.25 parabolic
6 0.5 -0.25 parabolic
Optimization terminated successfully:
the current x satisfies the termination criteria using
OPTIONS.TolX of 1.000000e-012
x =
0.5000
În cazul în care utilizăm subrutina fminsearch.m, iniţializarea
procesului de căutare se face printr-un singur punct şi nu printr-un interval
ca în cazul precedent. Pentru aceeaşi problemă ca şi în exemplu precedent
lansarea căutării se va face în forma:
xm=fminsearch(@f1,10,options)
unde 100 =x reprezintă condiţia iniţială de lansare a căutării. Rezultatele
sunt prezentate în cele ce urmează: Iteration Func-count min f(x) Procedure
1 2 90 initial
2 4 72 expand
3 6 42 expand
4 8 6 expand
5 10 2 reflect
6 12 0 contract inside
...
46 103 -0.25 shrink
47 106 -0.25 shrink
Optimization terminated successfully:
the current x satisfies the termination criteria using
OPTIONS.TolX of 1.000000e-012
and F(X) satisfies the convergence criteria using
OPTIONS.TolFun of 1.000000e-008
xm =
0.5000
9
În cazul funcţiilor de mai multe variabile, aplicaţiile se soluţionează
similar dar apar dificultăţi legate de înscrierea fişierelor funcţie.
Considerăm funcţia:
( )( )212121 sinsinsin),( xxxxxxf +++−= ,
pentru care se cere stabilirea unui minim pentru ]2/,0[, 21 π∈xx .
Fişierul funcţie fv3.m asociat funcţiei considerate este de forma:
function y=fv3(v)
x1=v(1);
x2=v(2);
y=-(sin(x1)+sin(x2)+sin(x1+x2));
Pentru a vizualiza graficul 3D al funcţiei, înscriem următorul fişier script:
clear all
map(1,:)=[rand(1) rand(1) rand(1)];
[x1,x2]=meshgrid(0:pi/60:pi/2);
z=-(sin(x1)+sin(x2)+sin(x1+x2));
surf(x1,x2,z)
xlabel('x1')
ylabel('x2')
Graficul funcţiei este prezentat în figura de mai jos:
00.5
11.5
2
00.5
1
1.52
-3
-2.5
-2
-1.5
-1
-0.5
0
x1x2
Figura 1. Graficul funcţiei
10
Dacă considerăm ca punct de iniţializare a căutării ( )0,00 =x , procesul de
căutare se lansează în forma:
xm=fminsearch(@fv3,[0,0],options)
Rezultatele sunt prezentate în continuare:
Iteration Func-count min f(x) Procedure
1 3 -0.0005 initial
2 5 -0.0015 expand
3 6 -0.0015 reflect
...
91 191 -2.59808 shrink
92 192 -2.59808 reflect
93 196 -2.59808 shrink
Optimization terminated successfully:
the current x satisfies the termination criteria using
OPTIONS.TolX of 1.000000e-012
and F(X) satisfies the convergence criteria using
OPTIONS.TolFun of 1.000000e-008
3. Chestiuni de studiat
Se consideră ca funcţie obiectiv una dintre funcţiile prezentate în tabelul de
mai jos (funcţia va fi impusă de conducătorul lucrării):
Nr. Crt. Funcţie Interval
1 22 012 1 15 1 55( ) . . .f x x x= ⋅ − ⋅ + [-3;3]
2 3 20 4375 2 012 2 437 1 552( ) . . . .f x x x x= − ⋅ + ⋅ + ⋅ + [-3;3] 3 4 3 20 4653 0 4375 6 465 2 438 1( ) . . . .f x x x x x= − ⋅ − ⋅ + ⋅ + ⋅ − [-3;3] 4 1 388 0 63980 368 2 235. .( ) . .x xf x e e− ⋅ ⋅= ⋅ + ⋅ [-2;2] 5 9 609 11 82 1 259 3 636 1 259( ) . . cos( . ) . sin( . )f x x x= − ⋅ ⋅ + ⋅ ⋅ [-2;2]
a) Determinaţi minimul funcţiei.
b) Afişaţi graficul funcţiei obiectiv.
11
Lucrare de laborator Metoda secţiunii de aur pentru minimizarea
funcţiilor de o singură variabilă
1. Scopul lucrării
Scopul principal al acestei lucrări este ca prin implementarea şi
rularea unui algoritm bazat pe metoda secţiunii de aur să familiarizeze pe cei
interesaţi cu problemele care apar în minimizarea funcţiilor de o singură
variabilă.
2. Prezentarea lucrării
Metoda “secţiunii de aur” constituie una dintre uneltele cele mai
puternice pentru minimizarea funcţiilor de o singură variabilă, bazaţi
exclusiv pe valorile funcţiei obiectiv.
Considerăm în cea de-a k-a iteraţie un interval de incertitudine [ak,bk].
Noul interval de incertitudine [ak+1,bk+1] se formează astfel:
⇒ [λk,bk] , dacă f(λk) > f(µk)
⇒ [ak, µk] , dacă f(λk) < f(µk)
În cadrul acestei metode, alegerea punctelor intermediare de lucru se
va face în următoarele impuneri:
i) Lungimea intervalului de incertitudine finală, în cadrul
oricărei iteraţii trebuie să fie independentă de rezultatul
comparaţiei f(λk) < f(µk) sau f(λk)> f(µk). Prin urmare:
bk – λ k = µk - ak (1)
2
12
O asemenea condiţie este indeplinită dacă alegem:
λ k = ak + (1-α )(bk-ak) (2)
cu )1,0(∈α şi
µk = ak + α ( bk-ak) (3)
Pentru o astfel de alegere, indiferent de rezultatul comparaţiei,
rezultă:
)1,0(,11 ∈=−− ++ αα
kk
kk
abab (4)
ii) Cea de-a doua condiţie impune ca în cadrul noii iteraţii, λ k+1 şi
µk+1 să rezulte λ k+1 = µk sau µk+1 = λ k. În aceste condiţii, în cea
de-a (k+1)-a iteraţie este necesară evaluarea unei singure
valori pentru funcţia obiectiv.
Pentru stabilirea completă a relaţiilor de definire (8) şi (9), este
necesar să stabilim valoarea α pentru care este îndeplinită cea de a doua
condiţie.
Presupunem că în cea de-a k-a iteraţie, rezultatul comparaţiei valorilor
funcţiei obiectiv este f(λk) > f(µk) şi, prin urmare, ak+1 = λ k , iar bk+1 = bk. În urma
condiţiei ii) impuse mai sus, este necesar ca µk+1 = λ k. Prin urmare:
λ k+1 = ak+1 + (1 - α )(bk+1 - ak+1) = λ k + (1-α )α (bk - ak) = µk . (5)
ak λk µk bk
ak+1 µk+1 bk+1
ak+1 λk+1 bk+1
λk+1
µk+1 dacă f(λk) > f(µk)
dacă f(λk) < f(µk)
Figura 2. Schemǎ de determinare a noului interval de cǎutare
13
Ţinând cont de modul de definire al valorilor λ k şi µk obţinem:
)()()1()()1( kkkkkkkk abaababa −⋅+=−⋅⋅−+−⋅−+ αααα (6)
şi în final, după prelucrări elementare:
012 =−+αα . (7)
Ecuaţia (7) pune în evidenţă singura soluţie interesantă:
)1,0(2
15∈
−=α .
În mod similar, se demonstrează că această valoare asigură
îndeplinirea condiţiei ii) şi pentru cazul f(λk) < f(µk).
Metoda propusă impune algoritmul “secţiunii de aur”:
Etapa de iniţializare. Se alege lungimea intervalului de incertitudine finală,
0>ε . Fie [a1,b1] intervalul de incertitudine iniţial, calculăm
λ 1 = a1 + (1-α )(b1-a1) şi µ1 = a1 + α ( b1 - a1),
în care
618.02
15≅
−=α .
Se calculează f(λ1) şi f(µ1), facem 1=k şi trecem la etapa de bază.
Etapa de bază. Pasul 1. Dacă bk –ak < ε , algoritmul se opreşte şi punctul
de optim este situat în intervalul [ak,bk]. Dacă nu, în ipoteza că f(λk) ≥ f(µk),
trecem la pasul 2, iar dacă f(λk) < f(µk) trecem la pasul 3.
Pasul 2. Considerăm ak+1 = λ k , bk+1 = bk , λk+1 = µk şi µk+1 = ak+1 +
α (bk+1 - ak+1). Calculăm valoarea funcţiei în f(µk+1) şi trecem la pasul 4.
Pasul 3. Considerăm ak+1 = a k , bk+1 = µk , µk+1 = µk şi λk+1 = ak+1 + (1 -
α )(bk+1 - ak+1). Calculăm valoarea funcţiei în f(λk+1) şi trecem la pasul 4.
Pasul 4. Se iterează 1+→ kk şi se reia de la pasul 1.
14
Evaluarea numărului de iteraţii necesare atingerii intervalului de
incertitudine finală se stabileşte cu uşurinţă. Pe baza relaţiei (6) stabilim:
n
nn
nn
abab
abab
abab
α=−−
⋅⋅−−
⋅−− ++ 11
22
33
11
22 (8)
prin urmare:
nnn
abab
α=−− ++
11
11 .
Condiţia de atingere a intervalului de incertitudine finală va fi:
εα <−⋅=− ++ )( 1111 abab nnn . (9)
Remarcăm că în cadrul acestei proceduri este necesar ca pe fiecare
iteraţie (cu excepţia primei iteraţii) se calculează o singură valoare a funcţiei
obiectiv (n+1 valori pentru asigurarea unei precizii ε cu n fixat de relaţia
(9)).
În continuare vom prezenta un exemplu privind utilizarea unui
asemenea algoritm pentru minimizarea unei funcţii. Funcţia obiectiv este de
forma:
1)( 2 ++= xxxf
şi admite un punct de minim pentru: 5.0min −=x şi 43)( min =xf .
Reamintim că în procedura de calcul a fost inserată o pauză de o
secundă pentru simularea unei complexităţi de calcul mai mare.
Pentru a evalua procedura de căutare propusă utilizăm subrutina
gold.m prezentată în Anexa 1. Subrutina a fost setată astfel:
• Intervalul de incertitudine iniţială 3,3 11 =−= ba ;
• Precizia de evaluare a minimului .01.0=ε
15
Rezultatele obţinute prin simulare sunt prezentate în tabelul următor:
k ka kb kk ab − )( kf λ )( kf µ
1 -3.0000 3.0000 6.0000 0.7933 2.2098
2 -3.0000 0.7082 3.7082 1.9242 0.7933
3 -1.5836 0.7082 2.2918 0.7933 0.8608
4 -1.5836 -0.1672 1.4164 1.0444 0.7933
5 -1.0426 -0.1672 0.8754 0.7933 0.7500
6 -0.7082 -0.1672 0.5410 0.7500 0.7659
7 -0.7082 -0.3738 0.3344 0.7565 0.7500
8 -0.5805 -0.3738 0.2067 0.7500 0.7522
9 -0.5805 -0.4528 0.1277 0.7510 0.7500
10 -0.5317 -0.4528 0.0789 0.7500 0.7503
11 -0.5317 -0.4829 0.0488 0.7502 0.7500
12 -0.5131 -0.4829 0.0301 0.7500 0.7500
13 -0.5131 -0.4944 0.0186 0.7500 0.7500
14 -0.5060 -0.4944 0.0115 0.7500 0.7500
15 -0.5016 -0.4944 0.0071 0.7500 0.7500
16 -0.5016 -0.4972 0.0044 0.7500 0.7500
17 -0.5016 -0.4988 0.0027 0.7500 0.7500
18 -0.5005 -0.4988 0.0017 0.7500 0.7500
19 -0.5005 -0.4995 0.0010 0.7500 0.7500
20 -0.5005 -0.4999 0.0006 0.7500 0.75
În prima coloană a acestui tabel este înscris tactul de lucru.
Următoarele două coloane indică valorile de capăt pentru intervalul de
incertitudine pe tactul curent. În coloana următoare sunt date lăţimea
intervalului de incertitudine pe tactul considerat. În sfârşit, ultimele două
coloane indică valorile funcţiei obiectiv în punctele intermediare de calcul
specifice metodei propuse.
Tabelul 2.1. Rezultatele obţinute prin simulare
16
Convergenţa algoritmului poate fi analizată conform graficului
prezentat în figura 3. Graficul prezintă evoluţia intervalului de incertitudine
în cadrul procesului de căutare a minimului. Astfel curba superioară indică
valorile maxime ale intervalelor de incertitudine iar cea inferioară indică
valorile minime a acestor intervale.
Lungimea intervalului de incertitudine pe tact este prezentată prin
segmentul trasat cu linia grasă.
Este interesantă o comparaţie a timpului de calcul necesar obţinerii
intervalului de incertitudine finală pentru diverse valori ale preciziei de
calcul impuse.
În tabelul 2.2 sunt prezentaţi timpii de calcul necesari stabilirii
intervalului de incertitudine finală pentru o precizie impusă 4,1,10 ∈− kk
pornind de la un interval de incertitudine iniţială .3,3 11 =−= ba
k 1.0000 2.0000 3.0000 4.0000
timp 11.1460 16.0230 21.0300 25.0360
Figura 3. Evoluţia intervalului de incertitudine
Tabelul 2.2. Variaţia timpului în funcţie de precizie
17
Graficul dependenţei timpului de calcul în funcţie de precizia impusă
este prezentat în figura următoare:
3. Chestiuni de studiat
Se consideră ca funcţie obiectiv una dintre funcţiile prezentate în tabelul de
mai jos (funcţia va fi impusă de conducătorul lucrării):
Nr. Crt. Funcţie
Interval de incertitudine
iniţial
Valoare de
minim
Minim
1 4 3 25 2( )f x x x x x= + + ⋅ + + [-3;3] -0.1027 1.9491
2 25 16( )f x x x= ⋅ + + [-4;4] -0.1 15.95
3 ( ) sin( ) cos( )f x x x= − [-2;2] -0.7854 -1.4142
4 5( ) cos( ) sin( )xf x e x x= − − ⋅ [-7;-3] -4.9112 -5.0917
5 2( ) cos( )f x x x x= − ⋅ [-3;4] 0.3889 -0.2086
Determinaţi minimul funcţiei, impunând diferite precizii de calcul al
acestuia. De asemenea, afişaţi graficul funcţiei obiectiv.
Figura 4. Evaluarea timpilor de calcul pentru diverse precizii impuse
18
Lucrare de laborator
Metoda căutării dihotomice
1. Scopul lucrării Scopul principal al acestei lucrări este ca prin implementarea şi
rularea unui algoritm bazat pe metoda căutării dihotomice să familiarizeze
pe cei interesaţi cu problemele care apar în minimizarea funcţiilor de o
singură variabilă.
2. Prezentarea lucrării
Procedurile de minimizare ce utilizează exclusiv valorile funcţiei, se
bazează pe o teoremă fundamentală în teoria funcţiilor convexe care permite
eliminarea din intervalul de incertitudine iniţială [a,b], subintervalul [µ,b]
sau [a,λ], în funcţie de rezultatul comparaţiei f(λ) şi f(µ). Dacă este
echiprobabilă eliminarea oricărui subinterval obţinem varianta optimă de
alegere a valorilor λ şi µ: 2
ab −== µλ , caz în care, în cadrul fiecărei
iteraţii eliminăm jumătate din intervalul de incertitudine iniţială. Din păcate,
o astfel de procedură nu poate fi aplicată datorită necesităţii comparării
valorilor f(λ) şi f(µ) care implică µλ ≠ . În cadrul metodei căutării
dihotomice, se impune alegerea valorilor λ şi µ în forma:
ηλ −−
=2
ab şi ηµ +−
=2
ab (1)
care apropie abordarea către varianta optimă în cazul η foarte mic.
3
19
În alegerea valorilor η trebuie să ţinem cont de zeroul maşinii care
limitează inferior ecartul de variaţie a acestei mărimi. De asemenea, trebuie
să ţinem cont de faptul că valorile funcţiei trebuie să fie comparabile, deci
valoarea de discernabilitate trebuie aleasă corespunzător tipului funcţiei şi al
maşinii utilizate.
Metoda căutării dihotomice generează următorul algoritm de căutare:
Etapa de iniţializare. Se alege constanta de discernabilitate η şi lungimea
intervalului de incertitudine finală 0>ε . Fie [a,b] intervalul de
incertitudine iniţială, facem 1=k şi trecem la etapa de bază.
Etapa de bază. Pasul 1. Dacă bk-ak < ε , algoritmul se opreşte, iar minimul
se găseşte undeva în intervalul [ak;bk]. Dacă nu, atunci se calculează:
ηλ −+
=2
kkk
ba , ηµ +
+=
2kk
k
ba şi se trece la pasul 2.
Pasul 2. Dacă )()( kk ff µλ < atunci kkkk baa µ== ++ 11 , iar în cazul în
care )()( kk ff µλ > , atunci kkkk bba == ++ 11 ,λ . Se iterează 1+→ kk şi se
reia pasul 1.
În privinţa numărului de iteraţii pentru o precizie finală impusă,
apreciem cu uşurinţă:
)211(2)(
21
1111 kkkk abab −⋅⋅+−⋅=− ++ η (2)
şi deci: εη <−⋅⋅+−⋅=− ++ )211(2)(
21
1111 nnnn abab
de unde putem determina numărul de iteraţii necesar obţinerii minimului cu
o precizie impusă.
20
Vom prezenta în continuare un exemplu simplu de minimizare a unei
funcţii de o singură variabilă utilizând metoda căutării dihotomice. Funcţia
ce urmează a fi minimizată este:
1)( 2 ++= xxxf
şi a fost înscrisă ca fişier cu extensia .m. Pentru a simula cazul unei funcţii
de calitate mult mai complexe şi care necesită un timp de calcul mai mare,
s-a forţat o pauză de 1 secundă care determină un timp relativ mare pentru
evaluarea funcţiei.
Soluţia de optim este evidentă .43,
21
=−= optopt yx Aplicăm
algoritmul căutării dihotomice pentru soluţionarea problemei enunţate
setând parametrii subrutinei prezentate în Anexa 1 astfel:
• Intervalul de incertitudine iniţială ]30,30[],[ 11 −=ba ; • Precizia de evaluare 01.0=ε ; • Valoarea de dicernabilitate 0001.0=e .
-3 -2 -1 -0.5 0 1 2 30
2
4
6
8
10
12
14
x
f(x)
x=-0.5punct de minim
Figura 5. Graficul funcţiei propuse ca exemplu
21
Valorile intermediare şi intervalul de incertitudine finală sunt
prezentate în tabelul urmǎtor:
k ka kb kk ab −
1.0000
2.0000
3.0000
4.0000
5.0000
6.0000
7.0000
8.0000
9.0000
10.0000
11.0000
12.0000
13.0000
14.0000
-30.0000
-30.0000
-15.0001
-7.5001
-3.7501
-1.8751
-0.9376
-0.9376
-0.7032
-0.5860
-0.5274
-0.5274
-0.5128
-0.5055
30.0000
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
-0.4686
-0.4686
-0.4686
-0.4686
-0.4979
-0.4979
-0.4979
60.0000
30.0001
15.0002
7.5002
3.7502
1.8752
0.9377
0.4689
0.2346
0.1174
0.0588
0.0295
0.0148
0.0075
Pentru a evidenţia corelaţia dintre timpul de calcul şi precizia metodei,
putem determina timpii necesari stabilirii intervalului de incertitudine finală
pentru diferite valori impuse preciziei de determinare a punctului de minim.
Astfel, impunând incertitudinea iniţială 5,5 11 =−= ba se determină intervalul
de timp necesar pentru a asigura o precizie de .4,3,2,1,10 ∈− kundek
k 1 2 3 4
timp 14.0200 20.0290 28.0400 34.90
Tabelul 3.1. Rezultatele metodei dihotomice
Tabelul 3.2. Variaţia timpului în funcţie de precizie
22
În tabelul 3.2 sunt prezentate rezultatele asupra testului propus, iar în
figura 6 este prezentată dependenţa grafică dintre timpul necesar procesării
şi precizia de calcul impusă.
Figura 6. Variaţia timpului de procesare în funcţie de precizie
3. Chestiuni de studiat
Se consideră ca funcţie obiectiv una dintre funcţiile prezentate în tabelul de
mai jos (funcţia va fi impusă de conducătorul lucrării):
Nr. Crt. Funcţie
Interval de incertitudine
iniţial
Valoare de
minim
Minim
1 22( ) sin( )f x x x x= ⋅ − ⋅ [-3;4] 0 0
2 27 31( )f x x x= ⋅ + + [-4;4] -0.0714 30.9643
3 2( ) cos( ) sin( )f x x x= − ⋅ [-1;4] 2.0345 -2.2361
4 3( ) cos( ) sin( )xf x e x x⋅= − − − [-4;1] -0.92 -1.3381
5 4 3 22 5 11( )f x x x x x= − ⋅ + ⋅ + + [-3;4] -0.0943 10.9519
Determinaţi minimul funcţiei, impunând diferite precizii de calcul al
acestuia. De asemenea, afişaţi graficul funcţiei obiectiv.
23
Lucrare de laborator
Metoda Fibonacci
1. Scopul lucrării
Scopul principal al acestei lucrări este ca prin implementarea şi
rularea unui algoritm bazat pe metoda Fibonacci să familiarizeze pe cei
interesaţi cu problemele care apar în minimizarea funcţiilor de o singură
variabilă.
2. Prezentarea lucrării
Metoda Fibonacci este o procedură de minimizare a funcţiilor de o
singură variabilă, în absenţa restricţiilor, bazată exclusiv pe valorile funcţiei
obiectiv. Ca şi metoda “secţiunii de aur”, metoda propusă caută minimizarea
numărului de valori calculate în cadrul procesului iterativ necesar asigurării
unei precizii impuse în localizarea minimului. Spre deosebire de metoda
“secţiunii de aur”, în cadrul acestei proceduri raportul de contracţie
kk
kk
abab
−− ++ 11
nu este constant şi este dependent de numărul iteraţiei efectuate. În
construcţia punctelor intermediare de test intervin valori ale termenilor din
şirul Fibonacci.
Şirul Fibonacci este definit prin relaţia de recurenţă: FK+1 = FK + FK-1 , cu F0 = F1 = 1. (1)
4
24
Dacă considerăm în cadrul iteraţiei k, intervalul de incertitudine
iniţială ca fiind [ak,bk], alegerea punctelor intermediare de test se face astfel:
)(1
1kk
kn
knkk ab
FF
a −⋅+=+−
−−λ (2)
)(1
kkkn
knkk ab
FF
a −⋅+=+−
−µ . (3)
În funcţie de rezultatul comparaţiei valorilor f(λk) şi f(µk), intervalul de
incertitudine va fi [ak, µk] sau [λk,bk]. Prin urmare:
)(1
11 kkkn
knkkkk ab
FF
aab −⋅=−=−+−
−++ µ (4)
sau
)()(11
111 kk
kn
knkk
kn
knkkkkkk ab
FF
abFF
abbab −⋅=−⋅−−=−=−+−
−
+−
−−++ λ . (5)
Prin urmare, indiferent de rezultatul comparaţiei valorilor f(λk) şi f(µk),
obţinem următorul rezultat:
)(1
11 kkkn
knkk ab
FF
ab −⋅=−+−
−++ . (6)
Urmează să verificăm dacă valorile intermediare calculate pe baza
relaţiilor (4.2) şi (4.3) satisfac condiţia pe baza căreia, la fiecare iteraţie, este
necesară calcularea valorii funcţiei într-un singur punct.
Considerăm că în iteraţia k, lungimea intervalului de incertitudine
iniţială este [ak,bk] şi alegem valorile λk < µk , conform relaţiilor (2) şi (3).
Presupunem cazul f(λk) > f(µk) care impune o incertitudine finală [ak+1 , bk+1]
= [λk , µk].
Aşadar:
25
=−⋅−−⋅+−⋅+=
=−⋅+=−⋅+=
+−
−−++
−
−−
+−
−−
−
−−++
−
−−++
))(()(
)()(
1
111
2
1
1
211
211
kkkn
knkk
kn
knkk
kn
knk
kkkn
knkkk
kn
knkk
abFF
abF
Fab
FF
a
abF
Fab
FF
a λλ
)()1()(1
12
1
1kk
kn
kn
kn
knkk
kn
knk ab
FF
FFab
FFa −⋅−⋅+−⋅+=
+−
−−
−
−−
+−
−− . (7)
Deoarece
11
11
1
11+−
−
+−
−−+−
+−
−− =−
=−kn
kn
kn
knkn
kn
kn
FF
FFF
FF
, (8)
rezultă:
kkkkn
knkkk
kn
kn
kn
knkk ab
FF
aabFF
FF
a µλ =−⋅+=−⋅⎟⎟⎠
⎞⎜⎜⎝
⎛++=
+−
−
+−
−−
+−
−−+ )()(
11
2
1
11 . (9)
Procedând prin calcul direct, obţinem pentru cazul f(λk) < f(µk):
kk λµ =+1 . (10)
Prin urmare, şi în cadrul acestei metode este necesară evaluarea
funcţiei obiectiv într-un singur punct, cu excepţia primei iteraţii.
Înainte de a prezenta algoritmic procedura de căutare, atragem
atenţia că iteraţiile se fac după valorile lui k pentru un n fixat apriori.
Valoarea n se determină pe baza preciziei impuse, ε. Dacă avem în vedere
relaţia (6), rezultă:
2
1
11 FF
abab
nn
nn =−−
−−
3
2
22
11
FF
abab
nn
nn =−−
−−
−− (11)
26
n
n
FF
abab 1
11
22 −=−−
Rezultatul se obţine imediat:
n
nn
Fabab 1
11
=−−
. (12)
Condiţia de asigurare a unei precizii impuse va fi:
)(111 ab
Fab
nnn −⋅=− (13)
din care putem determina valoarea lui n:
ε/)( 11 abFn −> .
În continuare, vom prezenta algoritmul Fibonacci:
Etapa de iniţializare. Se impune valoarea maximă a intervalului de
incertitudine finală ε, precum şi constanta de discernabilitate η. Se
consideră dat intervalul de incertitudine iniţială [a1,b1] şi determinăm n din
condiţia ε/)( 11 abFn −> . Se calculează:
)( 112
11 abF
Fa
n
n −⋅+= −λ , )( 111
11 abF
Fa
n
n −⋅+= −µ
şi valorile corespunzătoare: f(λ1) şi f(µ1).
Impunem 1=k şi trecem la etapa de bază.
Etapa de bază. Pasul 1. Dacă f(λk) > f(µk) se trece la pasul 2, iar dacă f(λk)
≤ f(µk) se trece la pasul 3.
27
Pasul 2. Se reiniţializează kka λ=+1 , kk bb =+1 şi se determină kk µλ =+1 ,
)( 112
11 ++−
−−++ −⋅+= kk
kn
knkk ab
FF
aµ .
Dacă 2−= nk , se trece la pasul 5 şi dacă nu, se trece la pasul 4.
Pasul 3. Se reiniţializează kk aa =+1 , kkb λ=+1 şi se determină:
)( 112
11 ++−
−−++ −⋅+= kk
kn
knkk ab
FF
aλ .
Dacă 2−= nk , se trece la pasul 5 şi dacă nu, se trece la pasul 4.
Pasul 4. Iterăm 1+→ kk şi reluăm de la pasul 1.
Pasul 5. Facem 1−= nn λλ şi ηλµ += nn . Dacă f(λk) > f(µk), atunci nna λ=
şi 1−= nn bb . În caz contrar, 1−= nn aa şi nnb λ= . Punctul de optim va fi
situat în intervalul [an , bn].
Exemplul 4.1. Comportarea algoritmului de căutare generat de metoda
Fibonacci va fi testatǎ pe funcţia 1)( 2 ++= xxxf , prezentată şi în
exemplele precedente şi pentru care am introdus o pauză de o secundă
simulând astfel o complexitate de calcul mai mare a funcţiei obiectiv.
Pentru căutarea minimului s-a utilizat algoritmul Fibonacci prezentat
în Anexa 1. Subrutina a fost setată astfel:
• Intervalul de incertitudine iniţială 3,3 11 =−= ba ;
• Precizia de determinare a minimului 01.0=ε ;
• Valoarea de discernabilitate 61 10−=e .
28
Rezultatele obţinute prin rularea programului sunt prezentate în tabelul 4.1:
k a(k) b(k) b(k)-
a(k) l(k) m(k) teta1(k) teta2(k)
1 -3.0000 3.0000 6.0000 -0.7082 0.7082 0.7933 2.2097
2 -3.0000 0.7082 3.7082 -1.5836 -0.7082 1.9242 0.7933
3 -1.5836 0.7082 2.2918 -0.7082 -0.1672 0.7933 0.8607
4 -1.5836 -0.1672 1.4164 -1.0426 -0.7082 1.0444 0.7933
5 -1.0426 -0.1672 0.8754 -0.7082 -0.5016 0.7933 0.7500
6 -0.7082 -0.1672 0.5410 -0.5016 -0.3738 0.7500 0.7659
7 -0.7082 -0.3738 0.3344 -0.5803 -0.5016 0.7565 0.7500
8 -0.5803 -0.3738 0.2066 -0.5016 -0.4525 0.7500 0.7523
9 -0.5803 -0.4525 0.1279 -0.5311 -0.5016 0.7510 0.7500
10 -0.5311 -0.4525 0.0787 -0.5016 -0.4820 0.7500 0.7503
11 -0.5311 -0.4820 0.0492 -0.5115 -0.5016 0.7501 0.7500
12 -0.5115 -0.4820 0.0295 -0.5016 -0.4918 0.7500 0.7501
13 -0.5115 -0.4918 0.0197 -0.5016 -0.5016 0.7500 0.7500
14 -0.5016 -0.4918 0.0098 -0.5016 -0.5016 0.7500 0.7500
Prima coloană indică tactul de lucru. Următoarele două coloane indică
limita inferioară şi respectiv superioară ale intervalului de incertitudine, iar
cea de a patra coloană indică lungimea intervalului de incertitudine pe tactul
curent. Coloanele cinci şi şase dau valorile de testare pe tactul curent, iar
ultimele două coloane fixează valorile funcţiei obiectiv în valorile
intermediare alese conform algoritmului propus.
Intuitiv, convergenţa algoritmului poate fi analizată cu ajutorul
graficului prezentat în figura 7.
Tabelul 4.1. Rezultatele obţinute prin simulare
29
Linia îngroşatǎ reprezintǎ intervalul de incertitudine pe fiecare tact.
Pentru o evaluare calitativă a vitezei de convergenţă a algoritmului
generat de metoda Fibonacci, pentru un acelaşi interval de incertitudine
iniţială am impus valori privind intervalul de incertitudine finală de forma
4,1,10 ∈= − kkε . Timpii de calcul necesari asigurării preciziei impuse sunt
prezentaţi în tabelul urmǎtor:
k 1.0000 2.0000 3.0000 4.0000
timp 12.2170 16.0430 21.0310 26.0370
În figura 8 este prezentat
graficul dependenţei timpului
de calcul în funcţie de precizia
de calcul impusă:
Figura 7.Convergenţa algoritmului
Tabelul 4.2. Variaţia timpului în funcţie de precizie
Figura 8. Dependenţa timpului în funcţie de precizia impusǎ
30
Observaţie
Algoritmul Fibonacci este superior ca performanţe algoritmului
“secţiunii de aur”. Acest lucru este resimţit în cazul unui număr mic de
iteraţii. În cazul unui număr mare de paşi, comportarea celor doi algoritmi
este aproape identică.
3. Chestiuni de studiat
Se consideră ca funcţie obiectiv una dintre funcţiile prezentate în tabelul de
mai jos (funcţia va fi impusă de conducătorul lucrării):
Nr. Crt. Funcţie
Interval de incertitudine
iniţial
Valoare de
minim
Minim
1 4 3 25 2( )f x x x x x= + + ⋅ + + [-3;3] -0.1027 1.9491
2 25 16( )f x x x= ⋅ + + [-4;4] -0.1 15.95
3 ( ) sin( ) cos( )f x x x= − [-2;2] -0.7854 -1.4142
4 5( ) cos( ) sin( )xf x e x x= − − ⋅ [-7;-3] -4.9112 -5.0917
5 2( ) cos( )f x x x x= − ⋅ [-3;4] 0.3889 -0.2086
Determinaţi minimul funcţiei, impunând diferite precizii de calcul al
acestuia. De asemenea, afişaţi graficul funcţiei obiectiv.
31
Lucrare de laborator
Căutarea liniară utilizând derivata funcţiei obiectiv
1. Scopul lucrării Scopul principal al acestei lucrări este ca prin implementarea şi rularea unui algoritm bazat pe metoda căutării liniare utilizând derivata funcţiei să familiarizeze pe cei interesaţi cu problemele care apar în minimizarea funcţiilor de o singură variabilă. 2. Prezentarea lucrării
Considerăm funcţia ℜ→],[:)( baxf şi acest interval [a,b] mărginit,
pe care funcţia este convexă şi deci, derivabilă. Dacă în cadrul iteraţiei k, iniţializată pe un interval de incertitudine [ak ,bk], alegem un punct de lucru
kkk ba << λ şi dacă evaluăm derivata )(| kf λ , atunci putem avea
următoarele situaţii:
i) 0)(| =kf λ . Prin urmare, kλ este un punct de staţionaritate şi
reprezintă soluţia de minim (deoarece )(xf este convexă);
ii) 0)(| >kf λ . În aceste condiţii, pentru kλλ >∀ , 0)()(| >−⋅ kkf λλλ
şi cum funcţia este convexă, rezultă: kkff λλλλ >∀≥ ),()( .
Aşadar, intervalul de incertitudine rezultat va fi în stânga lui kλ , deci:
[ak+1 , bk+1] = [ak , kλ ];
iii) 0)(| <kf λ şi 0)()(| >−⋅ kkf λλλ pentru kλλ <∀ şi prin urmare,
kkff λλλλ <∀≥ ),()( . Noul interval de incertitudine va fi în dreapta lui kλ :
[ak+1 , bk+1] = [ kλ , bk].
5
32
Aşadar, în funcţie de rezultatul comparaţiei, intervalul de incertitudine
după iteraţia k va fi [ak , kλ ] sau [ kλ , bk]. Alegerea valorii de test kλ se va
face astfel încât să minimizăm incertitudinea finală:
{ }),max(minarg kkkkoptimk ab −−= λλλ . (1)
Soluţia problemei de tip “minimax” este evidentă:
2
kkoptimk
ba +=λ . (2)
Valoarea de test se alege la jumătatea intervalului de incertitudine
[ak ,bk] şi indiferent de rezultatul comparaţiei,
2111 =
−− ++
kk
kk
abab
. (3)
Metoda impune imediat algoritmul cu înjumătăţirea intervalului de căutare:
Etapa de iniţializare. Fie [a1 , b1] intervalul de incertitudine iniţială şi ε
lungimea admisă pentru intervalul de incertitudine finală. Presupunem
cunoscute atât funcţia ℜ→],[:)( 11 baxf , cât şi derivata acesteia. Stabilim
cea mai mică valoare Ν∈n pentru care 11
121
abn −= , facem 1=k şi trecem
la etapa de bază.
Etapa de bază. Pasul 1. Considerăm 2
kkk
ba +=λ şi evaluăm )(| kf λ .
Dacă 0)(| =kf λ , algoritmul se opreşte şi kλ reprezintă soluţia de optim.
Dacă 0)(| >kf λ trecem la pasul 2, iar dacă 0)(| <kf λ trecem la pasul 3.
Pasul 2. Facem kk aa =+1 , kkb λ=+1 şi trecem la pasul 4.
33
Tabelul 5.1. Rezultatele obţinute prin simulare
Pasul 3. Facem kka λ=+1 , kk bb =+1 şi trecem la pasul 4.
Pasul 4. Iterăm 1+→ kk şi reluăm de la pasul 1.
În exemplul următor vom analiza concret comportarea algoritmului de înjumătăţire a intervalului de căutare utilizând derivata funcţiei obiectiv (în acest sens, vom considera cunoscută forma derivatei).
Exemplul 5.1. Funcţia obiectiv este 1)( 2 ++= xxxf , iar derivata se obţine
elementar în forma: 12 +⋅= xdxdf . În cadrul subrutinelor de evaluare a
funcţiei obiectiv şi a derivatei acesteia s-a inserat o pauză de o secundă pentru a simula o complexitate de calcul mai mare.
Pentru exemplificare s-a utilizat subrutina deriv1 setată astfel:
• Intervalul de incertitudine iniţială [ ]3,3 =−= ba ;
• Precizia de calcul 01.0=ε .
Rezultatele obţinute prin rularea programului sunt prezentate în tabelul 5.1:
k ka kb kk ab − 2/)( kk ba + ( )kcf ′ ( )kcf
1 -3.0000 3.0000 6.0000 0 1.0000 1.0000
2 -3.0000 0 3.0000 -1.5000 -2.0000 1.7500
3 -1.5000 0 1.5000 -0.7500 -0.5000 0.8125
4 -0.7500 0 0.7500 -0.3750 0.2500 0.7656
5 -0.7500 -0.3750 0.3750 -0.5625 -0.1250 0.7539
6 -0.5625 -0.3750 0.1875 -0.4688 0.0625 0.7510
7 -0.5625 -0.4688 0.0938 -0.5156 -0.0313 0.7502
8 -0.5156 -0.4688 0.0469 -0.4922 0.0156 0.7501
9 -0.5156 -0.4922 0.0234 -0.5039 -0.0078 0.7500
10 -0.5039 -0.4922 0.0117 -0.4980 0.0039 0.7500
11 -0.5039 -0.4980 0.0059 -0.5009 -0.0018 0.7500
34
Figura 9. Convergenţa algoritmului
Prima coloană indică tactul de lucru, iar următoarele două indică
capetele intervalului de incertitudine pe tactul curent. Cea de a patra coloană
indică lungimea intervalului lungimea intervalului de incertitudine pe tactul
curent iar cea de a cincea indică valoarea de înjumătăţire a tactului curent. În
sfârşit coloanele şase şi şapte fixează valorile pe care le ia derivata funcţiei
obiectiv şi funcţia obiectiv la jumătatea intervalului de incertitudine pe tact
curent.
Din tabel reiese că minimul se determină în 11 tacte, cu un interval de
incertitudine finală 498.0,5039.0[ 1111 −=−= ba ], eroarea de poziţionare
fiind egală cu .0059.0=ε Convergenţa către intervalul de incertitudine
finală poate fi urmărită în figura următoare:
Pentru o evaluare comparativă a timpului de calcul în dependenţă de
precizia de calcul impusă a fost rulată subrutina pentru acelaşi interval de
incertitudine iniţială dar precizii de calcul diferite, de forma
4,1,10 ∈= − kkε . Rezultatele obţinute sunt prezentate în tabelul 5.2:
35
Tabelul 5.2. Variaţia timpului în funcţie de precizie
k 1.0000 2.0000 3.0000 4.0000
timp 12.2480 16.0330 21.0400 26.0480
3. Chestiuni de studiat
Se consideră ca funcţie obiectiv una dintre funcţiile prezentate în tabelul de
mai jos (funcţia va fi impusă de conducătorul lucrării):
Nr. Crt. Funcţie
Interval de incertitudine
iniţial
Valoare de
minim
Minim
1 22 1( ) sin( )f x x x x= ⋅ − ⋅ + [-3;3] 0 1
2 27 31( )f x x x= ⋅ + + [-4;4] -0.0714 30.9643
3 2( ) cos( ) sin( )f x x x= − ⋅ [-1;4] 2.0345 -2.2361
4 3( ) cos( ) sin( )xf x e x x⋅= − − − [-4;1] -0.92 -1.3381
5 4 3 22 5 11( )f x x x x x= − ⋅ + ⋅ + + [-3;4] -0.0943 10.9519
6 2 2( ) sin( )f x x x x= − ⋅ + [-4;4] 0 2
7 1( ) ( cos( )) sin( )f x x x x= − ⋅ − − [-3;5] 3.4368 6.434
8 1( ) ( sin( )) cos( )f x x x x= ⋅ − + [1;5] -3 -4.4134
Determinaţi minimul funcţiei, impunând diferite precizii de calcul al
acestuia. De asemenea, afişaţi graficul funcţiei obiectiv.
36
Lucrare de laborator
Metoda căutării ciclice după axele de coordonate
1. Scopul lucrării Scopul principal al acestei lucrări este ca prin implementarea şi rularea unui algoritm bazat pe metoda căutării ciclice după axele de coordonate să familiarizeze pe cei interesaţi cu problemele care apar în minimizarea funcţiilor de o singură de mai multe variabile. 2. Prezentarea lucrării
Metoda pe care o prezentăm în continuare este o metodă rudimentară (metoda mai este cunoscută sub denumirea de “metodă naivă”), în care se exploatează succesiv posibilităţile de minimizare pe direcţia axelor de
coordonate nOxOxOx ,......., 21 , urmând ca procedura să fie reluată ciclic
până la atingerea (eventual) a unei zone de incertitudine finală admisă. În figura 3.1. este dată o reprezentarea geometrică a evoluţiei
procesului de căutare prin metoda propusă. În planul variabilelor 21Oxx ,
curbele trasate reprezintă aşa numitele curbe de izonivel care sunt grafice ale locului geometric din plan ce satisfac următoarea condiţie:
.),( 21 constxxf = (Din păcate, o astfel de reprezentare grafică este posibilă
numai pentru funcţii de două variabile).
Punctul de iniţializare a căutării se consideră ),( 21111 xxfy = .
Păstrând neschimbată valoarea 21x , urmărim o procedură de minimizare
după direcţia 1xO până la atingerea minimului parţial ),( 2112 xx . În
continuare, păstrăm valoarea 121 xx = procedând la o minimizare pe direcţia
6
37
2xO , până la atingerea valorii ),( 22212 xxfy = . În acest moment am încheiat
procedura de căutare pe un ciclu (am epuizat direcţiile oferite de axele de
coordonate) şi reluăm ciclul de căutare din punctul 2y .
Descrierea algoritmică a unei astfel de proceduri este următoarea:
Etapa de iniţializare. Se impune 0>ε condiţia de STOP (incertitudinea
finală). Se aleg direcţiile de căutare 11 ed = , 22 ed = , ... nn ed = , unde ei
reprezintă baza canonică ortogonală pentru nℜ . Se impune iniţializarea lui
1x şi se consideră 11 xy = şi 1== jk şi se trece la etapa de bază.
Etapa de bază. Pasul 1. Se determină
)(minarg1
jjj dyf ⋅+=ℜ∈
λλλ
şi se calculează jjjj dyy ⋅+=+ λ1 . Dacă nj < atunci 1+→ jj şi se reia
pasul 1. Dacă j = n, se trece la pasul 2.
Pasul 2. Considerăm 11 ++ = nk yx . Dacă ε<−+ |||| 1 kk xx , atunci algoritmul
se opreşte. În caz contrar, 11 += nxy , 1=j şi 1+→ kk şi se reia pasul 1.
Figura 10. Metoda căutării ciclice după axele de coordonate
),( 22212 xxfy =
),( 21111 xxfy = ),( 2112 xx
C1
C2<C1
C3<C2
C4<C3
x1
x2
0
38
Observaţii
Pentru o astfel de procedură, direcţiile de căutare sunt dictate de
modul de prezentare a funcţiei obiectiv. Oricare schimbare a axelor de
coordonate (rotaţii, translaţii) care nu schimbă efectiv minimul funcţiei şi
nici complexitatea, poate determina o creştere sau o scădere a eficienţei
metodei propuse (acest rezultat este exploatat în metoda căutării ciclice cu
pas accelerat ce urmează a fi prezentată).
Prezentăm în continuare modul de utilizare a acestei proceduri pentru
minimizare funcţiei:
( ) .3, 2132
3121 xxxxxxf ⋅⋅−+=
pornind de la un punct de iniţializare [ ] .7,70Tx = pentru o precizie .01.0=ε
-8 -6 -4 -2 0 2 4 6 8-8
-6
-4
-2
0
2
4
6
8
x1
x2
Figura 11. Curbele de izonivel pentru funcţia prezentată
39
Aplicând algoritmul prezentat materializat în subrutina coord.m,
obţinem următoarele rezultate: y =
1.0009
1.0005
A =
Columns 1 through 11
7.0000 2.6458 2.6458 1.2754 1.2754 1.0627
1.0627 1.0153 1.0153 1.0038 1.0038
7.0000 7.0000 1.6266 1.6266 1.1293 1.1293
1.0309 1.0309 1.0076 1.0076 1.0019
Columns 12 through 13
1.0009 1.0009
1.0019 1.0005
În figura 11 este prezentat graficul de evoluţie în procesul de căutare.
3. Chestiuni de studiat
Se consideră funcţia (forma funcţiei va fi impusă de conducătorul lucrării) a
cărei formă este prezentată în tabelul de mai jos:
Nr. Crt Funcţia Punct
de iniţializare
Punct deminim
Valoare de minim
1. 2 21 1 1 2 26 2 2 2x x x x x⋅ + ⋅ − ⋅ ⋅ + ⋅ [ ]1, 1− − [ ]2, 1− − -6
2. ( )221 2 1 2
1 109
x x x x+ + ⋅ + − [ ]1, 1− − [ ]5,0.5 7.5
3. ( ) ( )2 21 1 21 100x x x− + ⋅ − [ ]3, 4 [ ]1,1 0
4. ( ) ( )2 21 25 3 5x x⋅ − + − [ ]0,0 [ ]3,5 0
5. 2 21 1 2 2x x x x− ⋅ + [ ]1,2 [ ]0,0 0
40
6. 2 21 2 1 29 16 90 128x x x x⋅ + ⋅ − ⋅ − ⋅ [ ]0,3 [ ]5, 4 -481
7. 2 21 2 1 2 1 22 2 2 4 6+ + ⋅ − −x x x x x x [ ]1,1
1 4,3 3⎡ ⎤⎢ ⎥⎣ ⎦
143
−
8. 2 21 1 2 2 1 22x x x x x x− ⋅ + − ⋅ + [ ]3,5 [ ]1,0 -1
9. 21 1 2 2 1 25 4 16 12+ ⋅ + − −x x x x x x [ ]1,1 [ ]4,14− -152
10. 2 21 2 1 2 1 22 2 11 8+ + ⋅ − −x x x x x x [ ]3, 5− − [ ]2,3 -23
11. 2 21 2 1 1 2 22 2x x x x x x− + ⋅ + ⋅ ⋅ + [ ]1,1
31,2
⎡ ⎤−⎢ ⎥⎣ ⎦ -1.25
12. 2 21 2 1 2x x x x+ + ⋅ [ ]1,1 [ ]0,0 0
13. 2 21 216x x+ ⋅ [ ]2, 2 [ ]0,0 0
14. ( ) ( )2 21 1 21 x x x− + − [ ]5, 8− − [ ]1,1 0
a) Să se traseze graficul funcţiei considerate împreună cu curbele de
izonivel.
b) Să se evalueze eventualul punct de minim utilizând algoritmul propus.
c) Să se reprezinte grafic evoluţia în procesul de căutare.
41
Lucrare de laborator
Metoda căutării ciclice cu pas accelerat
1. Scopul lucrării Scopul principal al acestei lucrări este ca prin implementarea şi rularea unui algoritm bazat pe metoda căutării ciclice după axele de coordonate cu pas accelerat să familiarizeze pe cei interesaţi cu problemele care apar în minimizarea funcţiilor de o singură de mai multe variabile. 2. Prezentarea lucrării
Îmbunătăţirea procesului de căutare poate fi făcută intercalând între două etape de căutare ciclică după axele de coordonate, un proces de căutare
pe o direcţie impusă de soluţiile finale în cadrul ciclurilor amintite: kx şi
1+kx cu forma kkk xxd −= +1 .
În figura 12, pentru o funcţie obiectiv ipotetică, sunt prezentate primele două cicluri de căutare într-un algoritm cu
pas accelerat.
Căutarea se lansează în
1x şi după evaluările minimului în lungul direcţiilor 1Ox şi 2Ox , obţinem
punctul 2x . Urmează un “pas accelerat”; în fapt, un proces de minimizare pe
direcţia 121 xxd −= care oferă ca soluţie un punct y. Acest punct constituie
punctul de reiniţializare a unui nou procedeu de căutare ciclică după axele
x1
x2
x1
x2
x3
y y|
Figura 12. Un ciclu de cǎutare 0
7
42
de coordonate care stabileşte soluţia 3x . În continuare, intercalăm un nou
pas accelerat pe direcţia 232 xxd −= şi procedura continuă până la
atingerea condiţiei de STOP: ε<−+ |||| 1 kk xx .
Prezentăm în continuare algoritmic metoda propusă:
Etapa de iniţializare. Se impune 0>ε condiţia de încheiere a procesului
de căutare şi nx ℜ∈1 punctul de iniţializare a căutării. Considerăm 11 xy =
şi 1== jk şi se trece la etapa de bază.
Etapa de bază. Pasul 1. Calculăm
)(minarg1
jjj dyf ⋅+=ℜ∈
λλλ
şi construim jjjj dyy ⋅+=+ λ1 . Dacă nj < , atunci 1+→ jj şi se reia
pasul 1.
Dacă j = n, facem 11 ++ = nk yx . În condiţia ε<−+ |||| 1 kk xx , soluţia de optim
este 1+kx . În caz contrar, trecem la pasul 2.
Pasul 2. Construim direcţia de căutare kk xxd −= +1 şi determinăm
)(minarg 1*
1dyf k ⋅+= +
ℜ∈λλ
λ.
Impunem 1,1,*11 +→=⋅+= + kkjdxy k λ şi reluăm pasul 1.
Algoritmul prezentat este cunoscut în literatura de specialitate ca
“algoritmul Hooke şi Jeeves”. În multe probleme de minimizare a unor
funcţii convexe, utilizarea metodei căutării ciclice cu pas accelerat îşi
dovedeşte eficienţa oferind o aceeaşi precizie de determinare a minimului
într-un număr de paşi considerabil mai mic decât în cazul căutării ciclice
fără pas accelerat. Vom evidenţia acest avantaj pe următorul exemplu.
43
Considerăm funcţia obiectiv ce urmează a fi modificată de forma:
221
4121 )3()3(),( xxxxxf ⋅−+−= .
Pentru minimizarea funcţiei utilizăm metoda căutării ciclice cu pas
accelerat implementată prin subrutina Matlab sub denumirea coord_acc.
Pentru condiţia iniţială 10,10 21 == xx şi o precizie 001.0=e obţinem soluţia de minim:
0011.1,0032.3 21 == xx . Rezultatele comparative sunt prezentate în tabelul 7.1:
Cond
iniţiale Cătare ciclică Căutare ciclică
Pas
accelerat Căutare ciclică
10 5.3111 5.3111 3.8921 3.8921 3.0036 3.0032 3.0032
10 10 1.7704 1.7704 1.2973 1.0011 1.0011 1.0011
Atingerea soluţiei de minim se obţine în opt iteraţii în condiţiile
realizării preciziei de calcul impuse.
De remarcat este că în cazul unor condiţii iniţiale identice şi pentru o
aceeaşi precizie, convergenţa în punctul de minim se realizează în 113
iteraţii în cazul în care utilizăm metoda căutării ciclice după axele de
coordonate.
În figura 13 sunt prezentate evoluţiile procesului de căutare pentru
metoda căutării ciclice cu pas accelerat (vezi figura a) şi prin metoda
căutării ciclice după axele de coordonate (vezi figura b).
În continuare, vom analiza pe un exemplu concret modul în care
precizia de calcul impusă afectează timpii de calcul necesari asigurării unei
precizii impuse în cazul algoritmului de căutare ciclică cu pas accelerat.
Tabelul 7.1. Rezultatele pentru exemplul considerat
44
0 5 102
3
4
5
6
7
8
9
10
11
x1
x2
0 5 102
3
4
5
6
7
8
9
10
11
x1
x2
Pas accelerat
Considerăm funcţia (funcţia Rosenbrock):
21
221221 )1()(100),( xxxxxf −+−⋅= .
şi un punct de iniţializare 10,10 21 == xx .
Aplicăm procedura de minimizare utilizând metoda căutării ciclice cu
pas accelerat pentru o precizie de evaluare a punctului de minim de forma
3,2,1,10 ∈= − kkε şi vom stabili timpii de calcul necesari pentru a asigura
precizia impusă.
Precizia 1.0=ε 01.0=ε 001.0=ε
Căutare ciclică 0.07sec. 0.451sec. 47.067sec.
Căutare cu pas accelerat 0.12sec. 0.311sec. 4.36sec
a) cǎutare ciclicǎ cu pas accelerat b) cǎutare ciclicǎ dupǎ axele de coordonate
Figura 13. Rezultatele metodei dihotomice
Tabelul 7.2. Rezultatele pentru exemplu
45
Rezultatele sunt prezentate în tabelul 7.2, în care au fost incluse şi
datele referitoare la durata de calcul necesară asigurării aceleaşi precizii
utilizând metoda căutării ciclice după axele de coordonate.
1 2 30
10
20
30
40
50
Precizia impusa prin k.
Tim
pi d
e ca
lcul [s
ec.]
În figura de mai sus sunt prezentate dependenţele timpilor de calcul necesari
asigurării unei precizii impuse pentru cele două metode de căutare analizate.
Observaţii. O primă observaţie care trebuie făcută este legată de
posibilităţile de comparare a mai multor proceduri de căutare a minimului.
În principiu metodele pot fi comparate pe baza a trei criterii:
1. Precizia metodei evaluată prin distanţa dintre valoarea minimă a
funcţiei obiectiv obţinută prin aplicarea procedurii de căutare şi
valoarea de minim real sau prin distanţa dintre vectorul ce
caracterizează punctul de minim obţinut şi vectorul ce caracterizează
minimul real.
2. Numărul de evaluări numerice a funcţiei obiectiv (şi a derivatei
funcţiei obiectiv) necesar stabilirii punctului de minim.
3. Timpul de calcul necesar obţinerii soluţiei finale.
Figura 14. Dependenţa timpului de calcul - precizie
46
3. Chestiuni de studiat
Se consideră funcţia (forma funcţiei va fi impusă de conducătorul lucrării) a
cărei formă este prezentată în tabelul de mai jos:
Nr. Crt Funcţia
Punct de iniţializare
Punct deminim
Valoare de minim
1. 2 21 1 1 2 26 2 2 2x x x x x⋅ + ⋅ − ⋅ ⋅ + ⋅ [ ]1, 1− − [ ]2, 1− − -6
2. ( )221 2 1 2
1 109
x x x x+ + ⋅ + − [ ]1, 1− − [ ]5,0.5 7.5
3. ( ) ( )2 21 1 21 100x x x− + ⋅ − [ ]3, 4 [ ]1,1 0
4. ( ) ( )2 21 25 3 5x x⋅ − + − [ ]0,0 [ ]3,5 0
5. 2 21 1 2 2x x x x− ⋅ + [ ]1, 2 [ ]0,0 0
6. 2 21 2 1 29 16 90 128x x x x⋅ + ⋅ − ⋅ − ⋅ [ ]0,3 [ ]5,4 -481
7. 2 21 2 1 2 1 22 2 2 4 6+ + ⋅ − −x x x x x x [ ]1,1
1 4,3 3⎡ ⎤⎢ ⎥⎣ ⎦
143
−
8. 2 21 1 2 2 1 22x x x x x x− ⋅ + − ⋅ + [ ]3,5 [ ]1,0 -1
9. 21 1 2 2 1 25 4 16 12+ ⋅ + − −x x x x x x [ ]1,1 [ ]4,14− -152
10. 2 21 2 1 2 1 22 2 11 8+ + ⋅ − −x x x x x x [ ]3, 5− − [ ]2,3 -23
11. 2 21 2 1 1 2 22 2x x x x x x− + ⋅ + ⋅ ⋅ + [ ]1,1
31,2
⎡ ⎤−⎢ ⎥⎣ ⎦ -1.25
12. 2 21 2 1 2x x x x+ + ⋅ [ ]1,1 [ ]0,0 0
13. 2 21 216x x+ ⋅ [ ]2, 2 [ ]0,0 0
14. ( ) ( )2 21 1 21 x x x− + − [ ]5, 8− − [ ]1,1 0
a) Să se traseze graficul funcţiei considerate împreună cu curbele de izonivel.
b) Să se evalueze eventualul punct de minim utilizând algoritmul propus.
c) Să se reprezinte grafic evoluţia în procesul de căutare.
47
Lucrare de laborator
Algoritmul Rosenbrock
1. Scopul lucrării Scopul principal al acestei lucrări este ca prin implementarea şi
rularea unui algoritm bazat pe metoda Rosenbrock să familiarizeze pe cei
interesaţi cu problemele care apar în minimizarea funcţiilor de o singură de mai multe variabile.
2. Prezentarea lucrării
În cadrul acestui algoritm, direcţiile de căutare nddd …,, 21 care
constituie un sistem de vectori ortogonali, liniar independenţi se modifică
după fiecare ciclu de lucru. În general, procedura de lucru se lansează pe
direcţiile ( nddd …,, 21 ) fixate de axele de coordonate: Td )0,...,0,0,1(1 = , Td )0,...,0,1,0(2 = ,...., T
nd )1,0,...0,0(= .
Efectuând un ciclu de căutare complet pe aceste direcţii, obţinem:
)(minarg1
jjj dxf ⋅+=ℜ∈
λλλ
şi jjjj dxx ⋅+=+ λ1 .
După încheierea unui ciclu de căutare (toate direcţiile sunt
exploatate), se determină noi direcţii de căutare: **2
*1 ,...,, nddd astfel încât să
alcătuiască un sistem de vectori ortogonali, liniar independenţi şi care pentru
cazurile 0=iλ , sǎ îndeplineascǎ: ii dd =* .
Procedura de calcul este procedura de ortogonalizare Gram-Schmidt
aplicată astfel:
8
48
⎪⎩
⎪⎨
⎧
≠⋅
=
=∑=
n
jijii
jj
j dacad
dacada
0
0
λλ
λ ,
⎪⎩
⎪⎨
⎧
≥⋅⋅−
=
=∑−
=
1
1
** 2)(
1j
iii
Tjj
j
j jdacaddaa
jdacaab
||||
*
j
jj b
bd = (1)
Procedura de mai sus se poate scrie sub forma unui algoritm astfel:
Etapa de iniţializare. Se consideră dat 0>ε pentru evaluarea condiţiei de
STOP, nx ℜ∈1 punctul de iniţializare al căutării, sistemul de vectori
nddd …,, 21 (de obicei, axele de coordonate) care iniţiază procesul de
căutare. Cosiderăm 11 xy = şi 1== jk şi se trece la etapa de bază.
Etapa de bază. Pasul 1. Determinăm
)(minarg1
jjj dyf ⋅+=ℜ∈
λλλ
şi jjjj dyy ⋅+=+ λ1 .
Dacă nj < atunci 1+→ jj şi se reia pasul 1. În caz contrar, trecem la
pasul 2.
Pasul 2. Considerăm 11 ++ = nk yx . Dacă ε<−+ |||| 1 kk xx , algoritmul se
opreşte şi 1+kx este soluţia de optim.
Altfel, 11 += kxy , 1,1 +→= kkj şi trecem la pasul 3.
Pasul 3. Se construiesc conform (1) noile direcţii de căutare **2
*1 ,...,, nddd .
Facem 1*1 dd = , 2
*2 dd = ,..., nn dd =* şi reluăm de la pasul 1.
Algoritmul Rosenbrock, deşi mai complex în elaborare, are o eficienţă
sporită faţă de algoritmii anterior prezentaţi în condiţiile în care restricţiile
impuse pentru asigurarea convergenţei rămân aceleaşi.
49
Pentru prezentarea metodei Rosenbrock pentru câteva cazuri
concrete legate de minimizarea funcţiilor de mai multe variabile în absenţa
restricţiilor a fost implementat algoritmul în subrutina Matlab rosenbrock.
Se consideră funcţia RRxxf →221 :),( de forma:
( ) ( )221
4121 33),( xxxxxf ⋅−+−=
pentru care, considerând apriori că există un punct de minim (unic), se
impune determinarea acestuia utilizând algoritmul Rosenbrock.
Utilizând subrutina precizată anterior (rosenbrock.m) pentru, o
iniţializare arbitrară 2,3 21 == xx şi diverse valori ale lui ε (condiţia de
STOP) a fost testat algoritmul datele rezultate fiind prezentate în tabelul 8.1:
Precizia Soluţia de optim Valoare funcţiei obiectiv Număr de
iteraţii
1.0=ε 0001.10003.3
2
1
==
xx 12106748.3 −⋅ 9
01.0=ε 0001.10003.3
2
1
==
xx 12106748.3 −⋅ 9
001.0=ε 0001.10003.3
2
1
==
xx 12106748.3 −⋅ 9
0001.0=ε 0001.10002.3
2
1
==
xx 11107094.7 −⋅ 11
Se observă că pentru o precizie bună convergenţa este asigurată într-
un număr relativ mic de paşi iar valorile de optim obţinute prin rularea
programului sunt foarte apropiate de soluţia reală (aproape evident că soluţia
reală a problemei este 1,3 21 == xx cu un optim al funcţiei obiectiv
0=optf ).
Tabelul 8.1. Rezultatele pentru algoritmul Rosenbrock
50
Pentru condiţiile iniţiale specificate evoluţia în procesul de căutare
(pentru condiţia de stop 0001.0=ε ) este prezentată în figura 15, împreună
cu curbele de izonivel asociate funcţiei considerate:
-1 0 1 2 3 4 5 6-2
-1
0
1
2
3
4
0.1
1
10
100
100
10
200
20
0.1
110
100
100
1
200
20
x1
x2
Vom prezenta în continuare un exemplu de minimizare a unei funcţii
de trei variabile. Într-un astfel de caz nu mai există posibilitatea aprecierii
grafice a evoluţiei către punctul de minim (presupunând algoritmul
convergent); în schimb putem pune în evidenţă mult mai corect comportarea
algoritmilor în cazul unor probleme de mai mare complexitate.
Se consideră ),,( 321 xxxf : RR →3 de forma:
221
231
221321 )3()()2(),,( −++−+⋅−= xxxxxxxxxf .
Minimul funcţiei este asigurat de punctul de coordonate
2,1,2 321 === xxx la o valoare .0min =f
Figura 15. Curbele de izonivel şi evoluţia procesului de cǎutare
51
Ne propunem determinarea minimului funcţiei utilizând metoda
Rosenbrock. De asemenea se impune o comparaţie a acestei metode cu
celelalte metode de căutare a minimului fără ajutorul derivatelor studiate
până în acest moment:
Metoda Precizia Soluţia deoptim
Valoarea funcţieiobiectiv
Număr de iteraţii
1.0=ε 0279.20056.10279.2
3
2
1
===
xxx
0.0014 37
01.0=ε 0045.20009.10044.2
3
2
1
===
xxx
3.5612e-5 43
001.0=ε 0003.20001.10003.2
3
2
1
===
xxx
1.5051e-7 52
Cău
tare
cicl
ică.
0001.0=ε 0001.20000.10001.2
3
2
1
===
xxx
9.0768e-9 58
1.0=ε 9998.10000.19998.1
3
2
1
===
xxx
4.3982e-11 11
01.0=ε 9998.10000.19998.1
3
2
1
===
xxx
4.3981e-11 11
001.0=ε 9998.10000.19998.1
3
2
1
===
xxx
4.3981e-11 11
Cău
tare
cu
pas a
ccel
erat
0001.0=ε 0000.20000.10000.2
3
2
1
===
xxx
2.1523e-11 15
52
1.0=ε 0000.20000.10000.2
3
2
1
===
xxx
1.3081e-11 13
01.0=ε 0000.20000.10000.2
3
2
1
===
xxx
1.3081e-11 13
001.0=ε 0000.20000.10000.2
3
2
1
===
xxx
1.3081e-11 13
Ros
enbr
ock.
0001.0=ε0000.20000.10000.2
3
2
1
===
xxx
1.3081e-11 13
Pentru o analiză completă am considerat o iniţializare (altfel
arbitrară) 3,2,1.1000 == iptxi şi am rulat coord.m, coord_acc.m,
rosenbrock.m pentru diverse valori ale condiţiei de STOP.
Rezultatele obţinute sunt prezentate în tabelul 8.2.
Observaţii:
• Prezentarea rezultatelor trunchiate determină în unele cazuri
indicarea unor soluţii identice pentru precizii diferite. Într-o astfel de
situaţie observăm că valorile funcţiei obiectiv pentru o aceeaşi
soluţie (aparent) diferă. Este posibil însă ca impunând o precizie
relativ slabă soluţia obţinută să cadă într-un punct care asigură o
precizie mult mai bună decât cea impusă. În acest caz soluţiile pentru
două precizii distincte să fie identice iar valorile funcţiei obiectiv
sunt de asemenea identice.
Tabelul 8.2. Comparaţia celor 3 metode utilizate
53
• Metoda căutării ciclice după axele de coordonate este cea mai imprecisă şi totodată cea mai ineficientă.
• Deşi bazată pe o construcţie simplă metoda căutării ciclice cu pas accelerat se dovedeşte aproape la fel de precisă şi eficientă ca metoda Rosenbrock.
• Sub toate aspectele, metoda Rosenbrock este cea mai eficientă metodă din clasa celor analizate.
3. Chestiuni de studiat
Se consideră funcţia (forma funcţiei va fi impusă de conducătorul lucrării) a cărei formă este prezentată în tabelul de mai jos:
Nr.
Crt. Funcţia
Punct de
iniţializare
Punct
de
minim
Valoare
de
minim
1. 2 21 1 1 2 26 2 2 2x x x x x⋅ + ⋅ − ⋅ ⋅ + ⋅ [ ]1, 1− − [ ]2, 1− − -6
2. ( )221 2 1 2
1 109
x x x x+ + ⋅ + − [ ]1, 1− − [ ]5,0.5 7.5
3. ( ) ( )2 21 1 21 100x x x− + ⋅ − [ ]3, 4 [ ]1,1 0
4. ( ) ( )2 21 25 3 5x x⋅ − + − [ ]0,0 [ ]3,5 0
5. 2 21 1 2 2x x x x− ⋅ + [ ]1, 2 [ ]0,0 0
6. 2 21 2 1 29 16 90 128x x x x⋅ + ⋅ − ⋅ − ⋅ [ ]0,3 [ ]5, 4 -481
7. 2 21 2 1 2
1 2
2 2 24 6⋅ + ⋅ + ⋅ ⋅ −
− ⋅ − ⋅x x x x
x x [ ]1,1
1 4,3 3⎡ ⎤⎢ ⎥⎣ ⎦
143
−
54
8. 2 21 1 2 2 1 22x x x x x x− ⋅ + − ⋅ + [ ]3,5 [ ]1,0 -1
9. 2
1 1 2 2 1
2
5 4 1612⋅ + ⋅ ⋅ + − ⋅ −
− ⋅x x x x x
x [ ]1,1 [ ]4,14− -152
10. 2 21 2 1 2 1
2
2 2 118⋅ + ⋅ + ⋅ − ⋅ −
− ⋅x x x x x
x [ ]3, 5− − [ ]2,3 -23
11. 2 21 2 1 1 2 22 2x x x x x x− + ⋅ + ⋅ ⋅ + [ ]1,1
31,2
⎡ ⎤−⎢ ⎥⎣ ⎦ -1.25
12. 2 21 2 1 2x x x x+ + ⋅ [ ]1,1 [ ]0,0 0
13. 2 21 216x x+ ⋅ [ ]2,2 [ ]0,0 0
14. ( ) ( )2 21 1 21 x x x− + − [ ]5, 8− − [ ]1,1 0
a) Să se traseze graficul funcţiei considerate împreună cu curbele de izonivel.
b) Să se evalueze eventualul punct de minim utilizând algoritmul propus.
c) Să se reprezinte grafic evoluţia în procesul de căutare.
55
Lucrare de laborator
Tehnici de gradient
1. Scopul lucrării Scopul principal al acestei lucrări este ca prin implementarea şi rularea unui algoritm bazat pe metoda celei mai rapide coborâri să familiarizeze pe cei interesaţi cu problemele care apar în minimizarea funcţiilor de o singură de mai multe variabile.
2. Prezentarea lucrării
Considerăm funcţia obiectiv 1: ℜ→ℜnf , diferenţiabilă pentru nx ℜ∈∀ (f(x) este de clasa C1 pe ℜ ).
Definiţie. Fie un punct nkx ℜ∈ . Spunem că direcţia d constituie direcţie de
coborâre pentru f(x) în punctul kx dacă ( ) 0δ∃ > (oricât de mic) astfel că:
0( ) ( ), ( , )k kf x d f xλ λ δ+ ⋅ < ∀ ∈ .
În particular, definiţia prezentată, fixează condiţia “de coborâre” în forma:
0)()()(
lim00
<⋅∇=−⋅+
>→
dxfxfdxf
kTkk
λλ
λλ
. (1)
Atragem atenţia că modul de definire şi automat, condiţia de test (1),
au un caracter local deoarece sunt asociate unui punct de lucru kx .
Conform condiţiilor fixate privind funcţia obiectiv,
),()()()( xhOxfhxfhxf T +∇⋅=−+ , (2)
cu 0),(lim 1
0=⋅ −
→hxhO
h. Dacă 0)( ≠∇ xf atunci pentru h suficient de
mică, diferenţiala funcţiei f(x) poate fi scrisă în forma:
9
56
hxfxdf T ⋅∇= )()( . (3)
Conform inegalităţii Cauchy-Bunyakovski:
hxfhxfhxf T ⋅∇≤⋅∇≤⋅∇− )()()( . (4)
Pentru 0)( ≠∇ xf , partea dreaptă a inegalităţii (4) se transformă în
egalitate numai dacă )(xfh ∇⋅=α cu 0. >= constα . Pentru 0)( ≠∇ xf ,
partea stângă a inegalităţii se transformă în egalitate numai dacă )(xfh ∇⋅−= α cu 0. >= constα . Prin urmare, cea mai rapidă creştere se
obţine pentru evoluţia pe direcţia gradientului (într-o vecinătate a
lui nx ℜ∈ ). Pentru asigurarea celei mai rapide descreşteri într-o vecinătate a
lui nx ℜ∈ , evoluţia trebuie să se desfăşoare pe direcţia de antigradient.
Această proprietate a direcţiei de antigradient, deşi cum remarcam, are un caracter local, este deosebit de importantă şi stă la baza multor proceduri de minimizare a funcţiilor în absenţa restricţiilor. Una dintre aceste metode, cunoscută sub denumirea de metoda de gradient, va fi prezentată în cele ce
urmează. Se lansează problema pentru o iniţializare nx ℜ∈1 . Ca şi în cazul
algoritmilor anterior prezentaţi, nu există reţete pentru alegerea punctului iniţial de lucru. Dacă pe considerente ce ţin mai mult de natura procesului analizat, există o încadrare zonală pentru punctul de minim, este indicat ca iniţializarea să se facă cât mai aproape de punctul de minim.
Considerând că punctul de iniţializare a căutării a fost fixat, procedura
de gradient impune construcţia şirului }{ kx pe baza dependenţei recursive:
,...2,1,0,)(1 =>∇⋅−=+ kxfxx kkkkk λλ (5)
În relaţia (5), kλ reprezintă lungimea pasului (sau mai simplu, pasul)
metodei de gradient. Dacă 0)( ≠∇ xf , pasul kλ trebuie stabilit astfel încât
0)()( 1 <−+ kk xfxf . În baza relaţiei (2), putem evalua
57
0)()(
)()()()()( 1
<+⋅∇−=
=+∇⋅∇⋅−=−+
kkk
kkkT
kkk
OxfOxfxfxfxf
λλλλ
(6)
pentru kλ suficient de mic. Dacă 0)( =∇ xf , kx reprezintă un punct de
staţionaritate. În ipoteza că f(x) este o funcţie convexă, şi punctul de
staţionaritate este punct de minim. (În caz general, este necesară aprecierea
comportării funcţiei într-o vecinătate a lui kx , pentru a ne pronunţa dacă
acesta constituie sau nu un punct de minim).
Pentru alegerea lungimii pasului kλ sunt elaborate diverse metode,
generând la rândul lor diverşi algoritmi de căutare a minimului. Vom
prezenta câteva dintre principalele metode utilizate în practică.
1) Considerăm în cadrul iteraţiei k, punctul de iniţializare a căutării n
kx ℜ∈ şi )( kxf∇− direcţia de căutare (direcţia de anti gradient).
Considerăm funcţia de o variabilă α de forma:
0)),(()( ≥∇⋅−= ααλ kkk xfxfg (7)
Lungimea pasului de căutare se alege din condiţia:
0,))((inf)( *
0>=∇⋅−=
≥ kkkkkk gxfxfg ααλα
(8)
Metoda de gradient în care lungimea pasului de căutare se face
pe baza relaţiei (8) poartă numele de metoda celei mai rapide
coborâri.
Atragem atenţia asupra faptului că în general, stabilirea lui kα
pe baza relaţiei (8) nu este întotdeauna posibilă. Mai mult, este
posibil ca limita inferioară pentru (8) să nu fie accesibilă pentru
diverse valori ale indicelui k. În practică, kα se determină ca o
valoare care aproximează suficient de bine condiţia (8). În acest
sens, este posibil ca valoarea kα să se determine din condiţia:
58
∞<=≥+≤≤ ∑∞
=
δδδδα0
** ,0,)(k
kkkkkkk ggg (9)
sau
10,)0()1()( ** ≤<⋅+⋅−≤≤ kkkkkkkk gggg λλλα . (10)
Mărimile kk λδ , din (9) şi (10) caracterizează eroarea faţă de
condiţia (8); astfel încât, cu cât kδ se apropie mai mult de zero, iar
kλ de unitate, cu atât eroarea de evaluare în raport cu (8) este mai
mică.
Trebuie remarcat că direcţia de antigradient asigură o foarte
bună evoluţie doar în vecinătatea punctului de lansare a procedurii
de căutare. În aceste condiţii, dacă viteza de variaţie a funcţiei este
mare, atunci în punctul 1+kx , direcţia de antigradient ( 1( )kf x +−∇ )
poate diferi foarte mult de direcţia precedentă ( )( kxf∇− ). Din acest
motiv, stabilirea parametrului kα nu impune asigurarea relaţiei (8)
cu o mare fidelitate.
2) În practică, de multe ori ne limităm în alegerea pasului kα astfel
încât să asigurăm monotonia )()( 1 kk xfxf <+ fără a căuta
optimizarea pasului de căutare. Astfel, pentru iteraţia k:
)(1 kkkk xfxx ∇⋅−=+ λ . (11)
În cadrul fiecărei iteraţii, se evaluează modul în care evoluează
criteriul şi în funcţie de această evoluţie se modifică pasul astfel:
• dacă )()( 1 kk xfxf <+ , evoluţia este corespunzătoare şi se
încearcă modificarea în sens crescător în ideea accelerării
convergenţei (de exemplu, kk λλ ⋅=+ 5.11 );
59
• dacă )()( 1 kk xfxf =+ înseamnă că lungimea pasului de
căutare este prea mică încât diferenţele nu sunt sesizabile şi
ca şi în cazul precedent, mărim lungimea pasului de căutare;
• în sfârşit, dacă )()( 1 kk xfxf >+ , lungimea pasului de
căutare este prea mare. Aproximarea liniară care stă la baza
tehnicilor de gradient nu este corespunzătoare şi ca atare,
lungimea pasului trebuie diminuată ( kk λλ ⋅=+ 5.01 ) şi
rejectând valoarea obţinută 1+kx , reluăm iteraţia din kx .
3) Dacă funcţia obiectiv 1: ℜ→ℜnf este de clasă C1 şi dacă )( kxf∇
satisface condiţia Lipschitz:
( ) ( ) , , nf u f v L u v u v∇ −∇ < ⋅ − ∀ ∈ℜ , (12)
în care constanta L este cunoscută, atunci lungimea pasului de
căutare în tehnica de gradient poate fi aleasă de forma:
ε
αε⋅+
≤≤<220 0 Lk , (13)
în care εε ,0 reprezintă doi parametri de metodă. Dacă alegem
L/10 =ε şi 2/L=ε , obţinem metoda de gradient cu o lungime a
pasului de căutare ./1 constLk ==α Dacă constanta L este mare,
lungimea pasului va fi mică şi convergenţa (astfel asigurată) este
extrem de lentă.
Abordarea teoretică a problemei convergenţei algoritmilor de tip
gradient se porneşte de la premiza organizării procedurii iterative în forma:
0),(1 >∇⋅−=+ kkkkk xfxx αα , (14)
60
procedură nelimitată şi care conduce la un şir { kx }. În realitate, cum vom
vedea ulterior, algoritmul se opreşte în condiţia de STOP: ε<∇ )( kxf , cu ε
apriori impus. Considerând că succesiunea de valori formează un şir { kx },
problema care se pune este dacă acest şir este convergent către soluţia de
minim a problemei iniţiale, deci dacă se asigură convergenţa către mulţimea
punctelor de minim:
)}(inf)(,{ ** xffxfxXnx
n
ℜ∈==ℜ∈= (15)
sau altfel spus, dacă sunt îndeplinite condiţiile:
* *lim ( ) , lim ( , ) 0k kk kf x f x Xρ
→∞ →∞= = . (16)
Asigurarea convergenţei impune pe lângă condiţiile impuse prin
formularea problemei ( 1: ℜ→ℜnf , de clasă C1) o serie de restricţii
suplimentare destul de dure şi foarte greu de evaluat.
Fie ℜ→ℜ−∞>=ℜ∈
n
xxfxff
n:)(,)(inf* de clasă C1 şi pentru care
gradientul funcţiei satisface condiţia Lipschitz
nyxyxLyfxf ℜ∈∀−<∇−∇ ,,)()( .
În aceste condiţii, şirul { kx } construit iterativ în forma:
0,)(1 >∇⋅−=+ kkkkk xfxx αα , (17)
pentru care lungimea paşilor kα respectă:
∞<=≥+≤≤ ∑∞
=
δδδδα0
** ,0,)(k
kkkkkkk ggg , (18)
unde ))(()( kkk xfxfg ∇⋅−= αα pentru orice iniţializare nx ℜ∈1 , asigură
0)(lim =∇∞→ kk
xf .
61
În plus, dacă 1: ℜ→ℜnf este convexă şi dacă *f reprezintă
valoarea de minim atunci pentru { } )( 2−= kOkδ :
0.,)(0 01
0* >=⋅≤−≤ − constckcfxf k . (19)
Dacă considerăm că aceste restricţii sunt satisfăcute, condiţia de oprire
a algoritmului şi validarea soluţiei va fi ε<∇ )( kxf .
În continuare, este prezentat un algoritm de tip gradient în varianta
“celei mai rapide coborâri”:
Etapa de iniţializare. Fie 0>ε - condiţia de STOP. Alegem nx ℜ∈1
condiţia de iniţializare, 1=k şi trecem la etapa de bază.
Etapa de bază. Dacă ε<∇ )( kxf , algoritmul se opreşte şi kx reprezintă
soluţia de minim. Dacă nu, )( kk xfd −∇= şi se determină
))((minarg0
kkk xfxf ∇⋅−=>
λλλ
.
Construim )(1 kkkk xfxx ∇⋅+=+ λ , 1+→ kk şi se reia etapa de bază.
Interpretare geometricǎ. Putem da o interpretare geometrică intuitivă
legată de construcţia iterativă a soluţiei, utilizând algoritmul celei mai rapide
coborâri. Astfel, punctul 1+kx determinat în condiţiile:
)(1 kkkk xfxx ∇⋅−=+ λ şi ))((minarg0
kkk xfxf ∇⋅−=>
λλλ
se află pe dreapta { })(| kkkk xfxxxL ∇⋅−== λ în punctul de tangenţă a
acestei drepte cu conturul de izonivel
{ })()(| 11 ++ =ℜ∈=Γ kn
k xfxfx
şi dreapta kL este perpendiculară pe conturul de izonivel
{ })()(| kn
k xfxfx =ℜ∈=Γ .
62
Vom considera o parametrizare a curbei de izonivel kΓ de forma
)(txx = pentru bta ≤≤ şi deci .)())(( constxftxf k == pentru bta ≤≤ ,
iar kxtx =)( 0 . În aceste condiţii:
0))(()())(( =∇⋅= txftxtxfdtd T . (20)
În particular, pentru 0tt = obţinem:
0)()( 0 =∇⋅ kT xftx (21)
şi în felul acesta rezultă că direcţia de gradient (sau antigradient) este
perpendiculară pe kΓ în punctul kx . Reamintim că kλ se determină cu
condiţia:
0)()())(()()( 1 =∇⋅−∇=∇⋅−∇⋅−∇= +kkT
kkkT
kk xfxfxfxfxfg αα . (22)
Prin urmare, direcţia de gradient )( kxf∇ şi )( 1+∇ kxf sunt
perpendiculare şi cum )( 1+∇ kxf este perpendicular pe 1+Γk , rezultă că
)( kxf∇ este tangent curbei 1+Γk .
Considerentele geometrice prezentate sunt aplicate pe curbele de izonivel din figura 1.a şi 1.b:
Figura 16. Interpretarea geometrică a construcţiei iterative
63
Este evidentă o comportare mai bună pentru cazul 3.8.a.) în care
curbele de izonivel sunt apropiate unor contururi circulare. Pentru situaţia
prezentată în 16.b.), în care apar contururi alungite ce semnifică faptul că o
serie de variabile modifică mai puternic criteriul decât celelalte,
comportarea implică un zig-zag prelungit în care eficienţa căutării este
relativ mică.
Pentru a ilustra comportarea algoritmului celei mai rapide coborâri
vom prezenta în continuare câteva exemple.
Considerăm funcţia 221 :),( Rxxf :
221
4121 )3()3(),( xxxxxf ⋅−+−=
continuă şi derivabilă în 2R . Ne propunem stabilirea minimului acestei
funcţii utilizând tehnica de gradient. Evaluarea va fi făcută utilizând metoda
celei mai rapide coborâri implementată Matlab prin subrutina grad1.m.
Subrutina specificată a fost rulată pentru o iniţializare 10,10 21 == xx
pentru trei valori ale preciziei impuse .3,2,1,10 == − kkε Rezultatele
obţinute în urma rulării programelor precum şi durata procesării pentru
fiecare caz în parte sunt prezentate sintetic în tabelul 9.1:
Precizia Soluţia de minim
Valoarea funcţiei obiectiv
Timp necesarprocesării
Număr de iteraţii
110−=ε 0983.12937.3
2
1
==
xx 0074.0min =f 0.321 sec. 11
210−=ε 0442.11325.3
2
1
==
xx 4
min 1008.3 −⋅=f 0.861 sec.
27
310−=ε 0211.10634.3
2
1
==
xx 5
min 1061.1 −⋅=f 1.883 sec 81
410−=ε 0098.10295.3
2
1
==
xx 7
min 105433.7 −⋅=f 8.352 sec 335
Tabelul 9.1. Rezultatele utilizǎrii tehnicii de tip gradient
64
Pentru exemplul considerat se observă o bună comportare în privinţa
corectei soluţionări pentru o precizie convenabilă. Astfel, pentru condiţia de
stop 310−=ε obţinem soluţia 0211.1,0634.3 21 == xx faţă de soluţia reală
1,3 21 == xx . Creşterea exagerată a preciziei conduce la o îmbunătăţire
nesemnificativă a rezultatelor în dauna creşterii semnificative a timpului de
calcul. În acest sens, în figura 17 este prezentat graficul dependenţei
timpilor de calcul în funcţie de precizia impusă.
Din grafic reiese clar că o precizie exagerată conduce la un timp de
calcul extrem de mare.
1 2 3 40
1
2
3
4
5
6
7
8
9
k
timp
proc
esar
e (s
ec.)
Figura 17. Dependenţa timp de calcul - precizie
În figura 18 este prezentată dependenţa distanţei punctului de lucru
în a k-a iteraţie în raport cu punctul de minim real. Cum precizam şi în
prezentarea teoretică este evident că procedura aplicată are o eficienţă
ridicată în faza iniţială a algoritmului când distanţa în raport de minimul real
este relativ mare.
Rafinarea rezultatelor pentru asigurarea preciziei se face în paşi mici
care însă consumă timp de calcul.
65
0 10 20 30 40 50 60 70 80 900
2
4
6
8
10
12
Numar de iteratii
Distanta fata de minim
Figura 18. Dependenţa punct de lucru - optim
Deficienţa algoritmului anterior remarcată poate fi vizualizată şi în
graficul care reprezintă evoluţia procesului de căutare a punctului de minim.
Pentru cazul analizat, în condiţiile unei precizii 310−=ε , graficul evoluţiei
este prezentat în figura 19. Se observă cu uşurinţă că după circa patru iteraţii
căutarea se face într-un zig-zag prelungit extrem de ineficient.
3 4 5 6 7 8 9 101
2
3
4
5
6
7
8
9
10
110
100
500
1e+0031.5e+003
1
10
100
5001e+003
1.5e+003
x1
x2
Figura 19. Graficul evoluţiei către optim
66
Considerăm în continuare funcţia de două variabile:
133),( 22
321
3121 −⋅−+−⋅= xxxxxxf ,
pentru care gradientul se evaluează imediat în forma:
⎟⎟⎠
⎞⎜⎜⎝
⎛
⋅−⋅−⋅
=∇2
22
21
21 6319
),(xx
xxxf .
Curbele de izonivel pentru funcţia considerată sunt prezentate în
figura 20. Din alura acestor curbe rezultă că funcţia considerată nu este
unimodală pentru nRx∈ prezentând mai multe puncte de extrem. Ne
propunem stabilirea unui minim local utilizând tehnica de gradient.
Subrutina grad1.m implementează algoritmul de căutare propus.
-0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 11
1.2
1.4
1.6
1.8
2
2.2
2.4
2.6
2.8
3
-5
-5
-4
-3
-3
-2
-2
-
x1
x2
Solutia deminim local
Figura 20. Curbele de izonivel pentru funcţia consideratǎ
67
Considerând condiţia iniţială 8.2,8.0 21 == xx pentru o precizie
01.0=ε algoritmul este convergent în 4 iteraţii in punctul de minim local
0000.2,3333.0 21 == xx . Iniţializarea impusă a fixat punctul de plecare
într-o regiune a spaţiului pentru care soluţia obţinută constituie un atractor
stabil. O altă alegere a punctului de lansare a algoritmului putea să mă
conducă la o altă soluţie de minim local sau să genereze blocarea
algoritmului.
În următorul exemplu, considerăm drept funcţie obiectiv funcţia
Rosenbrock:
( ) ( )21
221221 1100),( xxxxxf −+−⋅=
pentru care gradientul este de forma:
( ) ( )( ) ⎟⎟
⎠
⎞⎜⎜⎝
⎛
−⋅−−−⋅⋅−
=∇ 212
12121
21 20012400
),(xx
xxxxxxf
Pentru o iniţializare arbitrară 2,2 21 == xx şi o precizie impusă 01.0=ε
convergenţa este asigurată în 3929 iteraţii. Soluţia aproximată pentru
problema propusă este 98321.0,99158.0 21 == xx . Valoarea de minim
obţinută pentru funcţia obiectiv este 5100973.7 −⋅=optf . Evoluţia în
procesul de căutare este prezentată în figura 21.
Legat de problema prezentată, putem face câteva observaţii de
ordin general:
• Condiţia de stop propusă în tehnicile de gradient este cu mult mai
relevantă decât în cazul metodelor de căutare directă. În cazul
problemei propuse, pentru o precizie de 1% obţinem o soluţie foarte
apropiată de soluţia reală.
68
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2
-0.5
0
0.5
1
1.5
2
2.5
3
x1
x2
Figura 21. Evoluţia în procesul de cǎutare
• Pentru obţinerea unei soluţii convenabile numărul de iteraţii este în general cu mult mai mic în raport cu metodele de căutare directă.
• În cazul funcţiei considerate evoluţia în procesul de căutare urmează un traseu curios în cadrul căreia există puncte care aparent sunt situate la o distanţă mai mare faţă de optim în raport cu punctul iniţial.
Figura 22. Succesiunea punctelor de test
69
Acest lucru apare mai evident dacă reluăm algoritmul propus pentru
o iniţializare 5,5 21 == xx . Obţinem o aproximare a soluţiei de
optim 0013.1,0007.1 21 == xx în 11320 de iteraţii la o valoare a
funcţiei obiectiv .103060.4 7−⋅=optf Succesiunea punctelor de test
este prezentată în figura 22.
3. Chestiuni de studiat
Se consideră funcţia (forma funcţiei va fi impusă de conducătorul lucrării) a
cărei formă este prezentată în tabelul de mai jos:
Nr.
Crt. Funcţia
Punct de
iniţializare
Punct
de
minim
Valoare
de
minim
1. 2 21 1 1 2 26 2 2 2x x x x x⋅ + ⋅ − ⋅ ⋅ + ⋅ [ ]1, 1− − [ ]2, 1− − -6
2. ( )221 2 1 2
1 109
x x x x+ + ⋅ + − [ ]1, 1− − [ ]5,0.5 7.5
3. ( ) ( )2 21 1 21 100x x x− + ⋅ − [ ]3, 4 [ ]1,1 0
4. ( ) ( )2 21 25 3 5x x⋅ − + − [ ]0,0 [ ]3,5 0
5. 2 21 1 2 2x x x x− ⋅ + [ ]1, 2 [ ]0,0 0
6. 2 21 2 1 29 16 90 128x x x x⋅ + ⋅ − ⋅ − ⋅ [ ]0,3 [ ]5, 4 -481
7. 2 21 2 1 2
1 2
2 2 24 6⋅ + ⋅ + ⋅ ⋅ −
− ⋅ − ⋅x x x x
x x [ ]1,1
1 4,3 3⎡ ⎤⎢ ⎥⎣ ⎦
143
−
8. 2 21 1 2 2 1 22x x x x x x− ⋅ + − ⋅ + [ ]3,5 [ ]1,0 -1
9. 2
1 1 2 2 1
2
5 4 1612⋅ + ⋅ ⋅ + − ⋅ −
− ⋅x x x x x
x [ ]1,1 [ ]4,14− -152
70
10. 2 21 2 1 2 1
2
2 2 118⋅ + ⋅ + ⋅ − ⋅ −
− ⋅x x x x x
x [ ]3, 5− − [ ]2,3 -23
11. 2 21 2 1 1 2 22 2x x x x x x− + ⋅ + ⋅ ⋅ + [ ]1,1
31,2
⎡ ⎤−⎢ ⎥⎣ ⎦ -1.25
12. 2 21 2 1 2x x x x+ + ⋅ [ ]1,1 [ ]0,0 0
13. 2 21 216x x+ ⋅ [ ]2, 2 [ ]0,0 0
14. ( ) ( )2 21 1 21 x x x− + − [ ]5, 8− − [ ]1,1 0
a) Să se traseze graficul funcţiei considerate împreună cu curbele de izonivel.
b) Să se evalueze eventualul punct de minim utilizând algoritmul propus.
c) Să se reprezinte grafic evoluţia în procesul de căutare.
71
Lucrare de laborator
Metoda Newton
1. Scopul lucrării Scopul principal al acestei lucrări este ca prin implementarea şi rularea unui algoritm bazat pe metoda Newton să familiarizeze pe cei interesaţi cu problemele care apar în minimizarea funcţiilor de mai multe variabile. 2. Prezentarea lucrării
În lucrările precedente am prezentat metode de ordinul întâi – cum mai sunt cunoscute metode de minimizare a funcţiilor de mai multe variabile care utilizează numai derivate parţiale de ordinul unu pentru funcţia obiectiv.
În cadrul acestor metode, pentru stabilirea direcţiei de coborâre am utilizat numai partea liniară a dezvoltării în serie Taylor. Dacă funcţia obiectiv este diferenţiabilă şi permite ca pe lângă vectorul gradient )(xf∇ să
construim relativ simplu hessianul, putem utiliza pentru minimizarea funcţiei metode de ordinul doi, care se bazează pe aproximarea cuadratică a funcţiei obiectiv. Cum aproximarea pătratică constituie o mai bună aproximare în raport cu cea liniară, este de aşteptat ca aceste metode să sporească eficienţa procedurilor de calcul.
2.1. Metoda Newton
Considerăm funcţia obiectiv f(x): nℜ → nℜ , de clasă Cn, convexă în nℜ . Fie x1∈ nℜ un punct de iniţializare a cǎutǎrii. Considerând construcţia
soluţiei optime printr-un proces iterativ, vom presupune că în cadrul iteraţiei
10
72
k, iniţializarea este impusă de xk ∈ nℜ . În acest punct, putem dezvolta
funcţia obiectiv în forma:
)(0)()(H)(!2
1)()()()( 2kkkkkkk xxxxxxxxxxfxfxf −+−⋅⋅−+−⋅∇+= Τ
Τ (1)
Considerăm aproximarea pătratică pentru variabila )(xgk astfel:
)()(H)(!2
1)()()()()( kkkkkk xxxxxxxxfxfxfxg −⋅⋅−+−⋅∇=−= ΤΤ
(2)
şi definim soluţia *kx dată de condiţia:
)(inf)(, **kk
xkk
nk xgxgx
nℜ∈=ℜ∈ . (3)
Condiţia de iniţializare a căutării în cea de a )1( +k -a iteraţie va fi
construită în forma:
)1,0(),( *1 ∈−+=+ kkkkkk xx ααλλ (4)
1) Pentru 1=kα , ,...1,0=∀k obţinem *1 kk xx =+ conform (1). Prin
urmare, iniţializarea iteraţiei (k + 1) implică soluţionarea (1). În cazul
funcţiilor convexe, condiţia (1) implică minimizarea aproximării cuadratice
şi deci:
0)()(H)()( 11 =−⋅+∇=∇ ++ kkkkkk xxxxfxg (5)
De fapt acest lucru impune ca pe fiecare iteraţie (3.29) cu 1=kα să
soluţionăm un sistem de ecuaţii algebrice (3) în raport cu ( kk xx −+1 ).
Dacă matricea H(xk) este nesingulară, atunci:
)()(11 kkkk xfxHxx ∇⋅−= −+ (6)
73
Dificultăţile de calcul sunt evidente. Trebuie ca pe fiecare iteraţie să
calculăm:
vectorul gradient )( kxf∇ ;
matricea hessian H(xk);
inversa matricei hessiane H-1(xk) (inversa unei matrice (n x n)
implică n3 înmulţiri).
Din acest motiv, această procedură se recomandă numai dacă vectorul
gradient şi matricea hessian pot fi calculate relativ simplu şi dacă problema
parţială de minimizare (1) poate fi simplu soluţionată.
Pe de altă parte, dacă metoda poate fi aplicată, viteza de convergenţă
este mult mai mare în raport cu metodele anterior prezentate.
2) Se impune în relaţia (2) kα în forma oik λα = , unde io reprezintă
cel mai mic număr pentru care este îndeplinită inegalitatea:
)())(()( **kk
ikk
ik xfxxxfxf ⋅⋅≥−+− ελλ (7)
unde 0>λ şi 1<ε sunt parametri de metodă.
3) Parametrul kα poate fi ales în forma:
10 ≤≤ kα )(min)(10
ααα kkk hh≤≤
= , cu ))(()( *kkkk xxxfh −+= αα (8)
Totuşi, dacă procedura este convergentă, rezultatele obţinute sunt
foarte bune şi sunt obţinute într-un număr mic de iteraţii ceea ce implică un
timp de procesare redus.
Pentru evaluarea comportării algoritmului Newton a fost creată
subrutina Matlab newton.m. Utilizând această subrutină, vom prezenta în
continuare următoarea aplicaţie.
74
Se consideră funcţia RRxxf →221 :),(
4 21 2 1 1 23 3( , ) ( ) ( )f x x x x x= − + − ⋅ ,
pentru care funcţia gradient este de forma:
⎟⎟⎠
⎞⎜⎜⎝
⎛
⋅−⋅−⋅−+−⋅
=∇)3(6
)3()3(4),(
21
213
121 xx
xxxxxf .
Este evident că funcţia considerată admite un minim unic în punctul
1,3 21 == xx pentru o valoare a funcţiei obiectiv .0min =f
Se propune analiza comportării algoritmului Newton direct pentru
minimizarea funcţiei considerate. Considerăm o iniţializare (altfel arbitrară)
10,10 21 == xx şi pentru diverse valori ale parametrului ε care impune
condiţia de stop rulăm subrutina newton.m. Rezultatele obţinute în urma
rulării programului sunt prezentate în tabelul 10.1. Deşi compararea
algoritmilor de căutare pe exemple concrete nu este concludentă, am
efectuat rulări în condiţii identice utilizând metoda celei mai rapide
coborâri:
Metoda Precizia Soluţia de optim
Valoarea funcţiei obiectiv
Număr de
iteraţii
1.0=ε 9598.08795.2
2
1
==
xx 2.1072e-4 8
01.0=ε 9598.08795.2
2
1
==
xx 2.1072e-4 8
001.0=ε 9827.09482.2
2
1
==
xx 7.1986e-6 79
Met
oda
New
tun
0001.0=ε 0024.10071.3
2
1
==
xx 2.5359-9 340
75
1.0=ε 0983.12937.3
2
1
==
xx 0.0074 11
01.0=ε 0442.11325.3
2
1
==
xx 3.0856e-4 27
001.0=ε 0211.10634.3
2
1
==
xx 1.6176e-5 81
Met
oda
de
grad
ient
0001.0=ε 0098.10295.3
2
1
==
xx 7.5433e-7 335
Într-o primă apreciere, observăm că algoritmul Newton funcţionează
corect, realizând convergenţa către o soluţie convenabilă într-un număr
relativ mic de iteraţii. Se observă o comportare mai bună a algoritmului
Newton în raport cu metoda celei mai rapide coborâri.
3. Chestiuni de studiat
Se consideră funcţia (forma funcţiei va fi impusă de conducătorul lucrării) a
cărei formă este prezentată în tabelul de mai jos:
Nr. Crt. Funcţia Punct de
iniţializarePunct deminim
Valoare de
minim
1. 2 21 1 1 2 26 2 2 2x x x x x⋅ + ⋅ − ⋅ ⋅ + ⋅ [ ]1, 1− − [ ]2, 1− − -6
2. ( )221 2 1 2
1 109
x x x x+ + ⋅ + − [ ]1, 1− − [ ]5,0.5 7.5
3. ( ) ( )2 21 1 21 100x x x− + ⋅ − [ ]3, 4 [ ]1,1 0
4. ( ) ( )2 21 25 3 5x x⋅ − + − [ ]0,0 [ ]3,5 0
Tabelul 10.1. Rezultatele obţinute pentru fiecare tip de metodǎ
76
5. 2 21 1 2 2x x x x− ⋅ + [ ]1, 2 [ ]0,0 0
6. 2 21 2 1 29 16 90 128x x x x⋅ + ⋅ − ⋅ − ⋅ [ ]0,3 [ ]5,4 -481
7. 2 21 2 1 2
1 2
2 2 24 6⋅ + ⋅ + ⋅ ⋅ −
− ⋅ − ⋅x x x x
x x [ ]1,1
1 4,3 3⎡ ⎤⎢ ⎥⎣ ⎦
143
−
8. 2 21 1 2 2 1 22x x x x x x− ⋅ + − ⋅ + [ ]3,5 [ ]1,0 -1
9. 2
1 1 2 2 1
2
5 4 1612⋅ + ⋅ ⋅ + − ⋅ −
− ⋅x x x x x
x [ ]1,1 [ ]4,14− -152
10. 2 21 2 1 2 1
2
2 2 118⋅ + ⋅ + ⋅ − ⋅ −
− ⋅x x x x x
x [ ]3, 5− − [ ]2,3 -23
11. 2 21 2 1 1 2 22 2x x x x x x− + ⋅ + ⋅ ⋅ + [ ]1,1
31,2
⎡ ⎤−⎢ ⎥⎣ ⎦ -1.25
12. 2 21 2 1 2x x x x+ + ⋅ [ ]1,1 [ ]0,0 0
13. 2 21 216x x+ ⋅ [ ]2,2 [ ]0,0 0
14. ( ) ( )2 21 1 21 x x x− + − [ ]5, 8− − [ ]1,1 0
a) Să se traseze graficul funcţiei considerate împreună cu curbele de izonivel.
b) Să se evalueze eventualul punct de minim utilizând algoritmul propus.
c) Să se reprezinte grafic evoluţia în procesul de căutare.
77
Lucrare de laborator Metoda gradienţilor conjugaţi.
Algoritmul Fletcher-Reeves
1. Scopul lucrării Scopul principal al acestei lucrări este ca prin implementarea şi rularea unui algoritm bazat pe metoda gradienţilor conjugaţi să familiarizeze pe cei interesaţi cu problemele care apar în minimizarea funcţiilor de mai multe variabile. 2. Prezentarea lucrării
Metoda gradienţilor conjugaţi este un caz particular al metodei generale a direcţiilor conjugate. Într-o primă formulare metoda a fost propusă de Hestenes şi Stiefel în 1952 şi utilizată pentru soluţionarea sistemelor algebrice liniare.
Vom considera pentru început că funcţia obiectiv este o funcţie pătratică, pentru care hessianul (care apare explicit în criteriu) este o matrice simetrică pozitiv definită.
În continuare vom prezenta o variantă a acestei metode concretizată în algoritmul Fletcher-Reeves (1964). În cadrul unei astfel de proceduri, direcţia de căutare pe iteraţia k este dată dintr-o combinaţie a vectorului antigradient şi a direcţiei de căutare în pasul anterior de forma:
11)( −− ⋅+−∇= kkkk dxfd β (1)
)( 11 xfd −∇=
Coeficientul 1−kβ se determină din condiţia ca dk să fie H-conjugat cu
vectorul dk-1. În aceste condiţii:
11
78
))((0 1111 −−Τ−
Τ− +−∇⋅=⋅⋅= kkkkkk dxfHddHd β . (2)
Rezultă cu uşurinţă:
11
11
)(
−Τ−
Τ−
− ∇⋅⋅∇⋅⋅
=kk
kkk dHd
xfHdβ . (3)
Următorul punct de lucru se obţine prin minimizare pe direcţia dk, pornind de la xk. Prin urmare:
kkkk dxx ⋅+=+ λ1 , (4)
unde )(minarg0
kkk dxf ⋅+=>
λλλ
.
Funcţia ce trebuie minimizată o considerăm de forma:
kkkkkkk ddxfdxfdxf ⋅⋅⋅+∇⋅⋅+=⋅+ ΤΤ H21)()()( 2λλλ . (5)
Minimul funcţiei de gradul doi, fixează lungimea pasului de căutare:
kk
kkk dd
xfd⋅Η⋅
∇⋅= Τ
Τ )(λ . (6)
Se poate demonstra cu uşurinţă că pentru o funcţie pătratică
xxxcxf Η+⋅= ΤΤ
21)( , vectorul gradient este de forma xcxf Η+=∇ )( .
Prin urmare:
kkkkkk
kkkkk
dxfdxcdxcxcxf
⋅Η⋅+∇=⋅Η⋅+⋅Η+==+Η+=⋅Η+=+∇ +
λλλ
)()()()1( 1 (7)
Din relaţia de definire a direcţiilor de căutare (1) rezultă:
)()()())(()( 11 kkkkkkkk xfxfxfdxfxfd ∇⋅−∇=∇⋅⋅+−∇=∇⋅ ΤΤ−−
ΤΤ β . (8)
Relaţia (7) a fost obţinută ţinând cont cǎ punctul xk se obţine prin
minimizarea funcţiei pe direcţia dk-1 şi prin urmare: 0)( =∇ kxf .
79
Pe de altă parte,
)())(( 1 kkkkkkkk xfHddxfddd ∇⋅⋅−=⋅+−∇⋅Η⋅=⋅Η⋅ Τ−−
ΤΤ β (9)
deoarece dk-1 şi dk sunt H-conjugate.
Prin urmare, dacă introducem (8) şi (9) în (5) obţinem:
)(
)()(
kk
kkk xfd
xfxf⋅Η⋅∇⋅∇
−= Τ
Τ
λ . (10)
În aceste condiţii, putem pune în evidenţă un rezultat interesant şi
anume că vectorii )( 1+∇ kxf şi )( kxf∇ sunt ortogonali.
Rezultǎ:
0)(
)()()(
)()(
))(()()()(
=⋅Η⋅∇⋅∇⋅Η⋅∇⋅∇
−∇⋅∇=
=⋅Η⋅+∇⋅∇=∇⋅∇
Τ
ΤΤ
ΤΤ
kkkk
kkkk
kkkkkk
dxfxfdxfxf
xfxf
dxfxfxfxf λ (11)
Elementele prezentate pun în evidenţă posibilitatea stabilirii a trei
şiruri de vectori: {dk} succesiunea direcţiilor de căutare, { )( kxf∇ }
succesiunea direcţiilor de gradienţi în punctele prin care evoluăm către
minim şi {xk} succesiunea punctelor regulate în procesul de căutare. Astfel,
pornind din condiţia de iniţializare x1 pe direcţia )( 11 xfd −∇= , determinăm
printr-un algoritm de optim unidimensional, x2. În continuare, evaluăm în
x2 vectorul gradient )( 22 xfd −∇= . Evident, )( 2xf∇ este ortogonal pe
direcţia d1.
În continuare, cunoscând )( 2xf∇ şi d1 fie baza relaţiei (1) se
determină d2, H-conjugat cu d1. vectorii )(),( 21 xfxf ∇∇ şi d1 sunt legaţi prin
(11). Apoi, procedura continuă în mod similar.
80
Apreciem imediat că succesiunea {dk} este un rezultat al procedurilor
de H-conjugare a vectorilor, aplicată gradienţilor { })( kxf∇ , iar succesiunea
{ })( kxf∇ se obţine prin ortogonalizarea vectorilor: ...,),( 211 QdQdxf∇
Mai departe, vom arăta că este necesar ca procedurile de obţinere a
vectorilor {dk}, { })( kxf∇ prin H-conjugarea sau ortogonalizarea vectorilor
gradient să fie simultane. Vom utiliza un procedeu inductiv [7]: vectorii
{dk} construiţi pe baza relaţiei (3.35) sunt H-conjugaţi, iar vectorii { })( kxf∇
sunt ortogonali. Pentru k = 1, evident cele două afirmaţii sunt adevărate.
Considerăm că cele două afirmaţii sunt adevărate până la unu şi deci
direcţiile {d1, d2, …, dk} sunt H-conjugate, iar vectorii gradient
{ })(),...,(),( 21 kxfxfxf ∇∇∇ sunt ortogonali. Vom demonstra că ele sunt
valabile şi pentru k + 1. Astfel:
).)()()(
))(()()()( 1
kikki
kkkiik
dxfxfxf
dxfxfxfxf
⋅Η⋅∇⋅+∇⋅∇=
=⋅Η⋅+∇⋅∇=∇⋅∇ΤΤ
Τ+
Τ
λ
λ (12)
Deoarece )( kxf∇ este ortogonal lui )( ixf∇ , pentru 1...,,2,1 −=∀ ki
şi cum conform (1): 0)()( 1 =∇⋅∇ +Τ
kk xfxf , obţinem:
kikik dxfxfxf ⋅Η∇⋅=∇⋅∇ Τ+
Τ )()()( 1 λ , (13)
pentru ki ,1∈∀ .
11)( −− ⋅+−=∇ iii ddxf β (14)
şi prin urmare:
kiiikik dddxfxf ⋅Η⋅⋅+⋅=∇⋅∇ Τ−−
Τ+
Τ )()()( 111 βλ . (15)
81
Direcţia dk este H-conjugată în raport cu di pentru oricare 1,1 −∈ ki .
Deoarece 0)()( 1 =∇⋅∇ +Τ
kk xfxf (vezi (1)) rezultă
0)()( 1 =∇⋅∇ +Τ
ik xfxf ki ,1∈∀ . (16)
Prima proprietate a fost astfel demonstrată. În privinţa direcţiilor de
căutare, considerăm că {d1, d2, …, dk} sunt H-conjugate şi va trebui să
arătăm că 01 =⋅Η⋅Τ+ ik dd pentru oricare ki ,1∈ . Pentru i = k, prin însăşi
modul de construcţie a direcţiilor de căutare 01 =⋅Η⋅Τ+ kk dd . Fie, pentru
11 +≤≤ ki :
)())(( 111 +Τ
+Τ
+Τ ∇⋅Η⋅−=⋅+−∇⋅Η⋅=⋅Η⋅ kikkkiki xfddxfddd β (17)
Ţinând cont de relaţia (16), obţinem: ))()(( 1 iiii xfxfd ∇−∇−=⋅Η⋅ +λ ,
în care 0≠iλ ( 0=iλ numai dacă 0)( =∇ ixf caz în care xi este punct de
staţionare). Prin urmare:
)())()((1111 +
Τ+
Τ+
Τ ∇⋅Η⋅∇−∇⋅−=⋅Η⋅ kiii
ki xfxfxfddλ
. (18)
Cum am demonstrat anterior, gradienţii sunt ortogonali şi deci:
01 =⋅Η⋅ +Τ
ki dd pentru ki ,1∈ . (19)
Algoritmul bazat pe metoda Fletcher-Reeves este prezentat în
continuare.
Etapa de iniţializare. Se impune 0>ε condiţia de STOP şi punctul de
iniţializare a căutării ni Rx ∈ . Considerăm y1 = x1, 1),( 11 ==−∇= jkxfd
şi trecem la etapa de bază.
82
Etapa de bază. Pasul 1. Dacă ε<∇ )( jyf algoritmul se opreşte şi yj
reprezintă optimul. În caz contrar, se determină
)(minarg0
jjj dyf λλλ
+=>
şi jjjj dyy λ+=+1 .
Dacă j < n se trece la pasul 2 şi dacă j = n trecem la pasul 3.
Pasul 2. Se stabileşte direcţia de căutare jjjj dyfd β+−∇=+ )(1 ,
unde 2
2
1
)(
)(
j
jj
yf
yf
∇
∇=
+β . Iterăm 1+→ jj şi reluăm pasul 1.
Pasul 3. Facem 1,1),(, 11111 +→=−∇=== ++ kkjyfdyxy kk şi
reluăm pasul 1.
În algoritmul prezentat intervin câteva modificări în raport cu metoda
expusă. Astfel:
i) În pasul 1, lungimea pasului se stabileşte prin căutare
unidimensională şi nu analitic. Această modalitate se impune pentru ca
metoda să poată fi extinsă şi în cazul în care f(x) este neliniară şi nu neapărat
pătratică.
ii) Stabilirea coeficientului jβ în pasul 2 apare modificată în raport cu
prezentarea făcută. Noua formă a coeficientului jβ este uşor deductibilă.
Construcţia matricelor de deflecţie
Am stabilit în paragraful precedent că stabilirea direcţiilor
H-conjugate {d1, d2, ..., dn} asociate unui criteriu de calitate pătratic:
λλλλ ⋅Η⋅⋅+⋅= ΤΤ
21)( cf , (20)
83
permite stabilirea punctului de minim în cel mult n – iteraţii (fiecare iteraţie
implică optimizarea unidimensională pe o direcţie impusă).
În continuare, va trebui să analizăm următoarele probleme:
i) stabilirea direcţiilor H-conjugate implică cunoaşterea matricei H
deci asocierea metodei strict către funcţii pătratice. O astfel de
restricţie este mult prea puternică şi este de dorit lărgirea cadrului în
care funcţionează aplicaţia;
ii) pentru cazul general, va trebui să stabilim o procedură de
construcţie a direcţiilor de căutare faţă a utiliza explicit matricea
hessian.
Legat de prima problemă facem observaţia că, în cazul unei iniţializări
în aproprierea punctului de minim, aproximarea pătratică este foarte bună.
În aceste condiţii, chiar dacă convergenţa nu va fi slabă în faza de
iniţializare, pe măsură ce mă apropii de minim aceasta creşte foarte mult
determinând pe ansamblu o funcţionare eficientă.
Ceea de a doua problemă ridicată, este de ordin constructiv şi va fi
soluţionată ca atare. Direcţiile de căutare se construiesc iterativ, în forma:
)( jj yfDd ∇⋅−= , (21)
prin urmare direcţiile de căutare vor fi direcţii de antigradient modificate
prin matricele de deflexie Dj convenabil alese astfel încât:
— direcţiile astfel construite să fie direcţii de coborâre;
— în cazul în care construcţia este aplicată pe funcţia obiectiv,
direcţiile obţinute în aplicarea algoritmului trebuie să fie H–conjugate. Acest
motiv determină definirea unor astfel de metode ca „metode quasi-newton”.
84
În cele ce urmează, vom prezenta într-o manieră simplificată tehnica
de construcţie a matricelor de defleţie.
Vom considera pentru început cazul unui criteriu de calitate pătratic,
prezentat în forma:
)()(21)()()()( kkkkk xxxxxxxfxfxf −⋅Η⋅−+−⋅∇+= ΤΤ , (22)
unde xk reprezintă punctul de iniţializare al căutării în iteraţia k, )( kxf∇ –
vectorul gradient în kλ . Matricea H (hessian) este simetrică şi pozitiv
definită (funcţia este pătratică şi deci H nu depinde de punctul de lucru).
Dacă xk+1 reprezintă punctul final obţinut în iteraţia k, vectorul gradient
evaluat în acest punct va fi:
)()()( 11 kkkk xxxfxf −⋅Η+∇=∇ ++ . (23)
Dacă iniţializarea căutării în iteraţia k se face în nR∈λ , propunem ca
pe această iteraţie matricea de deflecţie să fie Dk (aproximarea inversului
hessianului cu cadrul acestei iteraţii). Căutarea liniară impune.
)( 11 kkkkkk xxDxx −⋅Η=⋅−= ++ λ . (24)
Aproximarea inversei hessianului în noua iteraţie (matricea de
deflexie Dk+1) se construieşte în forma:
kkk CDD +=+1 , (25)
în care Ck este o matrice de corecţie simetrică de rangul său doi. Vom
prezenta modul de calcul a matricei de corecţie de rang unu şi vom relua
analiza şi pentru corectarea de rang doi care sunt mult mai eficiente.
Dacă Dk+1 reprezintă aproximarea universului hessianului, cum
kkkk xxxfxf −=∇−∇⋅ ++ 11-1 ))()((H ,
85
Dk+1 trebuie să satisfacă
kkkkk xxxfxfD −=∇−∇⋅ +++ 111 ))()((
şi deci
kkkkkk xxxfxfCD −=∇−∇⋅+ ++ 11 ))()(()( . (26)
Prin urmare
))()(()())()(( 111 kkkkkkkk xfxfDxxxfxfC ∇−∇−−=∇−∇ +++ . (27)
Dacă notăm direcţia de căutare în iteraţia k: )( kkk xfDd ∇⋅−= , din
(24) rezultă
kkkkk dxx µλ =⋅=−+1 . (28)
Pentru simplificare vom nota:
)()( 1 kkk xfxfq ∇−∇= + . (29)
Notaţiile introduse transpun relaţia (27) în forma
kkkkk
kkkkkkkkkkk q
qqDqDqD
qDqC ⋅⋅⋅−⋅−⋅−
=⋅−=⋅ Τ
Τ
)())((
µµµ
µ . (30)
Prin simplă identificare, din relaţia (30) rezultă termenul de corecţie în
forma:
kkkk
kkkkkkk qqD
qDqDD⋅−
−−= Τ
Τ
)())((
µµµ . (31)
Metoda se lansează cu D1=I. Prin urmare, prima iteraţie urmează o
banală tehnică a celei mai rapide coborâri. În cursul procesului iterativ,
matricele de deflexie aproximează tot mai bine inversa hessianului, iar
metoda se aproprie de performanţele metodei Newton.
86
În practică, sunt deosebit de mult utilizate corecţii prin matrici de rang
doi. Dintre aceste proceduri, foarte utilizate sunt metoda Davidon-Fletcher-
Powell (DFP), în care matricele de corecţie sunt generate în forma:
kkk
kkkk
kk
kkk qDq
DqqDq
C⋅⋅⋅⋅⋅
−⋅⋅
= Τ
Τ
Τ
Τ
µµµ (32)
şi metoda Broyden-Fletcher-Goldfarb-Shanna (BFGS), în care matricea de
corecţie este de forma:
kk
kkkkkk
kk
kk
kk
kkKK q
qDDqqq
qDqc Τ
ΤΤ
Τ
Τ
Τ
Τ ⋅⋅+⋅⋅−
⋅⋅
⋅⎟⎟⎠
⎞⎜⎜⎝
⎛⋅⋅⋅
+=µ
µµµ
µµµ
1 (33)
Observaţie. Comparativ, BFGS este mai sensibilă decât DFP în privinţa
erorilor care apar în optimizarea unidimensională pe fiecare pas de lucru.
Exemple
Prezentăm în cele ce urmează câteva exemple de utilizare a
algoritmului Fletcher-Reeves pentru minimizarea în absenţa restricţiilor a
unei funcţii de mai multe variabile. În acest scop, algoritmul generat pentru
metoda mai sus amintită a fost implementat ca funcţie Matlab sub
denumirea FReeves1.
Pentru cazul în care gradientul nu este cunoscut apriori, ci trebuie
evaluat prin program, a fost creată subrutina FReeves2, prezentată în aceeaşi
anexă.
Evaluarea algoritmului se face utilizând o aproximare de ordinul întâi.
Se consideră funcţia RRxxf →221 :),(
4 21 2 1 1 23 3( , ) ( ) ( ) ,f x x x x x= − + − ⋅
87
pentru care funcţia gradient este de forma:
⎟⎟⎠
⎞⎜⎜⎝
⎛
⋅−⋅−⋅−+−⋅
=∇)3(6
)3()3(4),(
21
213
121 xx
xxxxxf
Este evident că funcţia considerată admite un minim unic în punctul
1,3 21 == xx pentru o valoare a funcţiei obiectiv .0min =f
Se propune analiza comportării algoritmului Fletcher-Reeves pentru
minimizarea funcţiei considerate.
Alegem o iniţializare (altfel oricare) 15,100 21 == xx şi pentru o
precizie 001.0=ε lansăm subrutina FReeves1.m. Algoritmul este
convergent în 0467.31 =x , 0156.12 =x care aproximează suficient de bine
soluţia reală. Valoarea funcţiei obiectiv în punctul obţinut este
.107496.3 6min
−⋅≈f Convergenţa este asigurată în 32 de paşi.
Graficul evoluţiei în procesul de căutare este prezentat în figura 23:
Figura 23. Graficul evoluţiei în procesul de cǎutare
88
Pentru a ne forma o imagine mai clară în legătură de cu modul de
funcţionare a subrutinei am testat comportarea acesteia în condiţii identice
modificând numai condiţia de STOP (implicit precizia).
Rezultatele sunt prezentate sintetic în tabelul următor:
Precizia Soluţia de optim
Valoarea funcţieiobiectiv
Număr de iteraţii
1.0=ε 0911.12712.3
2
1
==
xx
0.0054 12
01.0=ε 0423.11267.3
2
1
==
xx 2.5747e-4
20
001.0=ε 0156.10467.3
2
1
==
xx
4.7496e-6 32
0001.0=ε 0095.10284.3
2
1
==
xx
6.5318e-7 40
Sub raport al timpului de procesare necesar elaborării soluţiei,
utilizând instrucţiunea tic-toc au fost obţinute următoarele rezultate:
Precizie 110−=ε 110−=ε 110−=ε 110−=ε
Timp de
procesare 5.658 sec 10.606 sec 18.346 sec. 23.894 sec.
Pentru a evidenţia dependenţa timpului de calcul de numărul de valori
ale funcţiei obiectiv ce trebuie evaluate în cursul procesǎrii, am crescut
artificial durata de calcul a funcţiei, introducând o pauzǎ de 10 milisecunde.
Tabelul 11.1. Rezultatele mai mulor rulǎri
Tabelul 11.2. Dependenţa timp de calcul - precizie
89
1 2 3 40
5
10
15
20
25
k (e=10e-k)
timp
(sec
.)
În figura 24 este prezentat graficul dependenţei timpului de calcul în
raport cu precizia impusǎ.
Pentru aceeaşi funcţie prezentată în exemplul precedent, să se testeze
comportarea algoritmului Fletcher-Reeves pentru cazul în care gradientul
este calculat on-line printr-o aproximare de ordinul întâi (subrutina
FReeves2.m).
Pentru o iniţializare 15,10 21 == xx (identică cu cazul prezentat în
exemplul precedent) s-au efectuat rulări ale subrutinei FReeves2.m pentru
diverse valori ale pasului utilizat în aproximarea gradientului şi pentru
diverse valori ale parametrului precizie.
Rezultatele rulării sunt prezentate în tabelul urmǎtor:
Pas pentru evaluarea
gradientului Precizie Soluţia de
optim Valoarea
funcţiei obiectivNumăr de
iteraţii
1.0=e 0605.11838.3
2
1
==
xx
0.0011 11 01.0=h
01.0=e Algoritmul nu funcţionează
Figura 24. Dependenţa timp - precizie
90
1.0=e 0913.12718.3
2
1
==
xx
0.0055 12
01.0=e 0361.11095.3
2
1
==
xx
1.4527e-4 15 001.0=h
001.0=e Algoritmul nu funcţionează
1.0=e 0932.12773.3
2
1
==
xx
0.0059e-4 12
01.0=e 0377.11129.3
2
1
==
xx
1.6245e-6 20 0001.0=h
001.0=e0051.10152.3
2
1
==
xx
5.2833e-8 28
Datorită faptului că gradientul este evaluat cu o aproximaţie
necontrolată şi deci, la rândul ei, oprirea algoritmului este de asemenea slab
controlată în rularea programului, apar rezultate aparent paradoxale. Astfel,
la un pas 01.0=h şi o precizie 1.0=ε , obţinem rezultate mai bune decât
pentru pasul 001.0=h şi o aceeaşi precizie. În orice caz, cel puţin pe
exemplul considerat, diferenţele între cele două proceduri sunt nesemni-
ficative. Vom vedea că, în privinţa timpului de calcul, diferenţele între cele
două proceduri sunt semnificative.
Precizie 1.0=ε 01.0=ε 001.0=ε
Evaluare gradientof-line
5.658 s 10.506 s 18.346 s.
Evaluare gradienton-line
8.533 s 15.131 s 22.503 s
Tabelul 11.3. Rezultate în urma mai multor rulǎri
Tabelul 11.4. Timpul de calcul necesar pentru diverse variante
91
În tabelul anterior au fost prezentaţi timpii de calcul necesari evaluării
soluţiei problemei formulate, în cele două conjuncturi:
utilizând evaluarea off-line a gradientului;
utilizând evaluarea on-line a gradientului.
Este evident că în cazul în care evaluăm gradientul on-line timpul de
calcul creşte considerabil.
Superioritatea algoritmului Fletcher-Reeves în raport cu algoritmii de
tip gradient apare în cazul unor funcţii de mai multe variabile de o mai mare
complexitate.
În cele ce urmează prezentăm o aplicaţie legată de minimizare unei
funcţii de trei variabile. Se consideră funcţia RRf →3: de forma:
221
231
221321 )3()()2(),,( −++−+⋅−= xxxxxxxxxf
având funcţia gradient de forma:
⎟⎟⎟
⎠
⎞
⎜⎜⎜
⎝
⎛
⋅+⋅−−⋅+⋅−⋅−⋅−⋅
=∇
31
21
321
321
226102
226),,(
xxxx
xxxxxxf .
Funcţia prezintă un minim unic
2,1,2 321 === xxx
pentru o valoare a funcţiei obiectiv prezentate: .0min =f
Problema propusă este de a determina minimul acestei funcţii
utilizând algoritmul Fletcher-Reeves (implementat Matlab) şi de a realiza o
comparaţie cu utilizarea algoritmului de tip gradient în condiţii identice.
92
În acest sens, am considerat condiţia de iniţializare a căutării
10,10,10 321 === xxx şi au fost rulate subrutinele FReeves1.m şi grad1.m
pentru diverse valori ale parametrului precizie.
Rezultatele obţinute pun în evidenţă, cel puţin pentru exemplu
considerat, superioritatea metodei Fletcher-Reeves asupra tehnicii celei mai
rapide coborâri atât în privinţa preciziei în stabilirea soluţiei cât şi a
numărului redus de iteraţii necesare stabilirii soluţiei:
Metoda Precizie Soluţia de minim
Valoare funcţiei obiectiv
Număr de iteraţii
1.0=ε 9998.10030.19998.1
3
2
1
===
xxx
6.7001e-5 4
01.0=ε 9974.19999.09993.1
3
2
1
===
xxx
4.6722e-6 8
Met
oda
Flet
cher
-Ree
ves
001.0=ε0000.20000.10000.2
3
2
1
===
xxx
9.0904e-12 10
1.0=ε 0695.29999.00348.2
3
2
1
===
xxx
0.0036 26
01.0=ε 0092.20000.10046.2
3
2
1
===
xxx
6.3079e-5 36
Met
oda
de g
radi
ent
001.0=ε0008.20000.10004.2
3
2
1
===
xxx
4.8612e-7 48
Tabelul 11.5. Rezultate. Comparaţie
93
Deşi mai puţin sugestivă faţă de cazul funcţiilor de două variabile, în
figura de mai jos este prezentat graficul evoluţiei procesului de căutare
pentru iniţializarea considerată şi precizia 0 01.ε = :
02
46
810
0
5
100
2
4
6
8
10
x1x2
x3
Solutia deminim.
Figura 25. Graficul evoluţiei procesului de cǎutare
Pentru evaluarea timpilor de procesare în ambele metode, am
modificat forma funcţiilor introducând în cadrul funcţiei un timp mort fictiv,
simulând astfel o complexitate sporită pentru funcţia obiectiv. Utilizând
instrucţiunea tic-toc, am evidenţiat timpii de calcul pentru obţinerea soluţiei
în ambele proceduri. Rezultatele sunt prezentate în tabelul urmǎtor:
Precizia 1.0=ε 01.0=ε 001.0=ε 0001.0=ε
Metoda
Newton 0.58 s 0.581 s 6.319 s 27.23 s
Metoda de
gradient 4.52 s 16.73 s 64.86 s 323.57 s.
Tabelul 11.6. Timpul de execuţie pentru fiecare metodǎ în parte
94
3. Chestiuni de studiat
Se consideră funcţia (forma funcţiei va fi impusă de conducătorul lucrării) a
cărei formă este prezentată în tabelul de mai jos:
Nr. Crt. Funcţia Punct de
iniţializarePunct de minim
Valoare de
minim
1. 2 21 1 1 2 26 2 2 2x x x x x⋅ + ⋅ − ⋅ ⋅ + ⋅ [ ]1, 1− − [ ]2, 1− − -6
2. ( )221 2 1 2
1 109
+ + ⋅ + −x x x x [ ]1, 1− − [ ]5,0.5 7.5
3. ( ) ( )2 21 1 21 100x x x− + ⋅ − [ ]3, 4 [ ]1,1 0
4. ( ) ( )2 21 25 3 5x x⋅ − + − [ ]0,0 [ ]3,5 0
5. 2 21 1 2 2x x x x− ⋅ + [ ]1, 2 [ ]0,0 0
6. 2 21 2 1 29 16 90 128x x x x⋅ + ⋅ − ⋅ − ⋅ [ ]0,3 [ ]5,4 -481
7. 2 21 2 1 2 1 22 2 2 4 6⋅ + ⋅ + ⋅ ⋅ − ⋅ − ⋅x x x x x x [ ]1,1
1 4,3 3⎡ ⎤⎢ ⎥⎣ ⎦
143
−
8. 2 21 1 2 2 1 22x x x x x x− ⋅ + − ⋅ + [ ]3,5 [ ]1,0 -1
9. 21 1 2 2 1 25 4 16 12⋅ + ⋅ ⋅ + − ⋅ − ⋅x x x x x x [ ]1,1 [ ]4,14− -152
10. 2 21 2 1 2 1 22 2 11 8⋅ + ⋅ + ⋅ − ⋅ − ⋅x x x x x x [ ]3, 5− − [ ]2,3 -23
11. 2 21 2 1 1 2 22 2x x x x x x− + ⋅ + ⋅ ⋅ + [ ]1,1 [ ]1,3 / 2− -1.25
12. 2 21 2 1 2x x x x+ + ⋅ [ ]1,1 [ ]0,0 0
13. 2 21 216x x+ ⋅ [ ]2,2 [ ]0,0 0
14. ( ) ( )2 21 1 21 x x x− + − [ ]5, 8− − [ ]1,1 0
a) Să se traseze graficul funcţiei considerate împreună cu curbele de izonivel.
b) Să se evalueze eventualul punct de minim utilizând algoritmul propus.
c) Să se reprezinte grafic evoluţia în procesul de căutare.
95
Lucrare de laborator
Metoda Davidon-Fletcher-Powell
1. Scopul lucrării Scopul principal al acestei lucrări este ca prin implementarea şi rularea unui algoritm bazat pe metoda Davidon, Fletcher şi Powell şi să continue analiza problemelor care pot apărea în minimizarea funcţiilor de mai multe variabile.
2. Prezentarea lucrării Metoda pe care o vom prezenta în cele ce urmează a fost introdusă
iniţial de Davidon (1959) şi dezvoltată de Fletcher şi Powell (1963). Pe considerente ce vor fi prezentate ulterior, metoda este cunoscută şi ca „metoda cu metrică variabilă”. Ea face parte din metodele de căutare
quasi-newton, metode în care direcţiile de căutare se evaluează în forma
)( jkj yfDd ∇−= . (1)
Direcţia de antigradient este modificată (rotită) prin multiplicare cu o
matrice nxnjD ℜ∈ simetrică şi pozitiv definită care aproximează pe fiecare
pas matricea hessian. Matricele Dj se construiesc iterativ astfel ca matricea Dj+1 se se obţinǎ prin însumarea matricei Dj cu doi termeni de corecţie formulaţi de matrice de rangul doi. Vom începe prin prezentarea algoritmului generat printr-o astfel de procedură.
Etapa de iniţializare. Fie 0>ε condiţia de stop. Se impune punctul de
iniţializare nRx ∈1 şi matricea D1 simetrică şi pozitiv definită (de obicei se
alege D1 = In). Considerăm y1 = x1, k = j = 1 şi se trece la etapa de bază.
12
96
Etapa de bază. Pasul 1. Dacă ε<∇ )( jyf , atunci algoritmul se opreşte şi
yj constituie soluţia. Dacă nu, construim direcţia de căutare
)( jjj yfDd ∇⋅−= şi jjjj dyy λ+=+1 , unde )(minarg0
jjj dyf λλλ
+=>
. Dacă
j < n se trece la pasul 2. Dacă j = n, facem 1,1,111 =+→== ++ jkkyxy nn
şi se reia pasul 1.
Pasul 2. Se construieşte matricea de deflexie Dj+1:
jjj
jjjj
jj
jjjj qDq
DqqDq
DD⋅⋅
⋅⋅⋅−
⋅
⋅+= Τ
Τ
Τ
Τ
+ µµµ
1 ,
unde:
).()( 1 jjj
jjj
yfyfq
d
∇−∇=
⋅=
+
λµ
Se reiniţializează 1+→ jj şi se reia pasul 1.
Am atras atenţia asupra faptului că în condiţiile în care ne aflăm la o
distanţă relativ mare în raport cu minimul, algoritmul nu are o comportare
foarte bună. În orice caz, într-o astfel de locaţie, direcţiile astfel calculate
trebuie să asigure o direcţie de coborâre.
În cele ce urmează, vom demonstra că matricele de deflexie D1, D2,
…, Dn sunt simetrice şi pozitiv definite cu implicaţia imediată că, d1, d2, …,
dn sunt direcţii de coborâre.
Astfel, prin condiţiile de generare, D1 este o matrice simetrică şi
pozitiv definită. Direcţia )( 111 yfDd ∇⋅−= este o direcţie de coborâre
deoarece:
0)()()( 11111 <∇⋅⋅−∇=⋅∇ ΤΤ yfDxfdyf . (2)
97
Presupunem, că pentru un jDDDnj ...,,1 ,2,1−≤ sunt simetrice şi
pozitiv definite şi vom arăta că Dj+1 are aceleaşi proprietăţi. Simetria nu
reprezintă o problemă şi ea poate fi verificată direct. Considerăm, pentru
0≠∀λ :
jjj
jjjj
jj
jjjj qDq
DqqDxq
xxxDxxDx
⋅⋅
⋅⋅⋅⋅−
⋅
⋅⋅⋅+⋅⋅=⋅⋅ Τ
ΤΤ
Τ
ΤΤΤ
+Τ
µµµ
1 . (3)
Matricea Dj este simetrică şi pozitiv definită şi prin urmare există
matricea 2/1jD , astfel ca 2/12/1
jjj DDD ⋅= . În aceste condiţii, dacă notăm
jjj qDbxDa ⋅=⋅= 2/12/1 , expresia (3) devine:
jj
jj q
pxbb
babbaaxDx⋅
⋅+
⋅⋅−⋅⋅
=⋅⋅ Τ
Τ
Τ
ΤΤΤ
+ µ
22
1
)()())(( . (4)
Utilizǎm inegalitatea Cauchy-Bunyakovski-Schwartz:
0)()()( 2 ≥−⋅ ΤΤΤ babbaa
şi pentru că 01 ≥+Τ xDx j , este necesar ca: bTb > 0 şi 0>Τ
jj qp .
Cu notaţiile introduse:
))()(( 1 jjjjjj yfyfdq ∇−∇⋅=⋅ +ΤΤ λµ . (5)
Pasul de căutare jλ se alege prin minimizare unidimensională pe
direcţia dj şi deci 0)( 1 =∇⋅ +Τ
jj yfd . Aşadar,
)()()( jjjjjjjjj yfDyfyfdq ∇⋅⋅∇⋅=∇⋅⋅−=⋅ ΤΤ λλµ . (6)
Cum Dj este pozitiv definită şi 0)( ≠∇ jyf , rezultă cǎ 0>⋅Τjj qµ
( 0>jλ deoarece căutarea se face în lungul direcţiei de coborâre cu 0>λ ).
98
Rezultă că 0>jq şi deci, 0>⋅⋅=⋅Τjjj qDqbb . Prin urmare, am
demonstrat că 01 ≥⋅⋅ +Τ xDx j . Considerăm 01 =⋅⋅ +
Τ xDx j . O asemenea
condiţionare impune ca simultan 0)()()( 2 =−⋅ ΤΤΤ babbaa şi 0=⋅Τ xjµ .
Reamintim că 0)()()( 2 =−⋅ ΤΤΤ babbaa dacă şi numai dacă vectorii a şi b
sunt coliniari, deci dacă ba λ= .
Condiţia ba ⋅= λ implică jjj qDxD ⋅⋅=⋅ 2/12/1 λ . Aşadar: jqx ⋅= λ .
Cum 0≠x , rezultă 0≠λ . Cea de-a doua cerinţă este 0=⋅Τ xjµ şi deci
0=⋅⋅ Τjj qµλ şi cum 0≠λ ar fi necesar ca 0=⋅Τ
jj qµ , infirmând
rezultatele obţinute anterior. Prin urmare, 0,01 ≠∀>⋅⋅ +Τ xxDx j şi deci
Dj+1 este pozitiv definită. Direcţia dj+1 pentru 0)( ≠∇ jyf se construieşte în
forma )( 111 +++ ∇−= jjj yfDd :
0)()()( 11111 <∇⋅⋅−∇=⋅∇ +++Τ
++Τ
jjjjj yfDyfdyf . (7)
Se poate demonstra că dacă algoritmul este aplicat unei funcţii pătratice, direcţiile construite pe principiul enunţat sunt H-conjugate şi deci
*1 λ=+ny fiind punctul de optim (şi în aceste condiţii, evident, 1
1 H−+ =nD ).
Pentru a putea evalua comportarea algoritmului Davidon-Fletcher-Powell au fost elaborate două subrutine de calcul bazate pe algoritmul menţionat. O primă subrutină implementează algoritmul în condiţiile în care gradientul funcţiei a fost calculat off-line şi deci putem dispune de forma analitică a acestuia. Această subrutină este denumită dfp.m şi este prezentată în Anexa 1. O a doua subrutină implementează acelaşi algoritm, în ipoteza în care gradientul funcţiei nu este apriori cunoscut şi acesta trebuie evaluat în cadrul programului (derivatele parţiale sunt calculate prin aproximări de ordinul unu utilizând subrutina nabla.m). Subrutina intitulată dfp1.m se aflǎ tot în anexa mai sus menţionată.
99
Se consideră funcţia:
( ) ( )221
4121 33),( xxxxxf ⋅−+−= ,
având gradientul uşor de precizat în forma:
⎟⎟⎠
⎞⎜⎜⎝
⎛
⋅−⋅−⋅−⋅+−⋅
=∇)3(6
)3(2)3(4),(
21
213
121 xx
xxxxxf .
Punctul de minim pentru funcţia considerată este evident
1,3 21 == xx care asigură valoarea de minim .0min =f Problema propusă
este de a evalua comportarea algoritmului Davidon-Fletcher-Powell în
stabilirea aproximativă a punctului de minim.
Tabelul 12.1. Rezultatele obţinute
S-a considerat un punct de iniţializare 5,4 21 == xx şi s-a rulat
subrutina dfp.m pentru diverse valori ale preciziei ε impusă condiţie de
stop. Rezultatele obţinute sunt prezentate sintetic în tabelul 12.1.
În figura 26 este prezentată evoluţia procesului de căutare pentru
iniţializare anterior precizată şi o condiţie de stop .001.0=ε
Precizia Soluţia de minim
Valoarea funcţiei obiectiv
Număr de iteraţii
1.0=ε 0723.12152.3
2
1
==
xx
0021.0min =f 6
01.0=ε 0402.10609.3
2
1
==
xx 4
min 100942.2 −⋅=f 8
001.0=ε 0203.10609.3
2
1
==
xx 5
min 103767.1 −⋅=f 14
0001.0=ε 0402.10289.3
2
1
==
xx 6
min 100163.7 −⋅=f 44
100
-1 0 1 2 3 4 5 6 7-3
-2
-1
0
1
2
3
4
5
0.1
1
10
100
200
200
300
300
300
3
0.1
1
10
100
200
200300
300
3
3
x1
x2
Figura 26. Evoluţia procesului de cǎutare
În tabelul 12.2 sunt prezentate valorile timpilor de procesare pentru
asigurarea preciziei impuse.
Precizia 1.0=ε 01.0=ε 001.0=ε 0001.0=ε
Timp de
procesare 0.120 sec. 0.441 sec. 0.441 sec 0.901 sec
În cel de al doilea exemplu, se consideră funcţia:
( ) ( )221
4121 33),( xxxxxf ⋅−+−= ,
având gradientul uşor de calculat în forma:
⎟⎟⎠
⎞⎜⎜⎝
⎛
⋅−⋅−⋅−⋅+−⋅
=∇)3(6
)3(2)3(4),(
21
213
121 xx
xxxxxf .
Tabelul 12.2. Timpul de procesare în funcţie de precizie
101
Funcţia precizată a fost minimizată în exemplul anterior utilizând
algoritmul Davidon-Fletcher-Powell considerând cunoscută forma
gradientului. Vom relua aceeaşi problemă de minimizare aplicând acelaşi
algoritm în condiţiile în care presupunem că gradientul este construit on-line
printr-o aproximare de ordinul unu (subrutina dfp1.m completată cu
nabla.m). Considerăm o aceeaşi iniţializare 3,4 21 == xx iar pasul utilizat
pentru aproximarea derivatelor parţiale ce compun gradientul se alege
.0001.0=h Pentru diferite valori ale parametrului ε care impune condiţia
de stop şi implicit precizia metodei rezultatele rulării sunt prezentate în
tabelul următor:
Precizia Soluţia de optim
Valoarea funcţiei obiectiv
Număr de iteraţii
1.0=ε 0696.12071.3
2
1
==
xx
0018.0min =f 6
01.0=ε 0166.10497.3
2
1
==
xx
6min 100871.6 −⋅=f 8
001.0=ε 0165.10497.3
2
1
==
xx
6min 100815.6 −⋅=f 9
0001.0=ε Algoritmul nu funcţionează
Cum precizam anterior algoritmul se comportă impredictibil, dar pe exemplul considerat şi pentru alegerea făcută pentru pasul de aproximare
algoritmul se comportă bine.
În privinţa timpului de calcul, cum era de aşteptat acesta este mai
mare decât în cazul în care funcţia gradient este cunoscută apriori. În tabelul
4 sunt prezentate valorile timpilor de calcul pentru diverse precizii de
evaluare a soluţiei.
Tabelul 12.3. Rezultatele obţinute
102
Precizia 1.0=ε 01.0=ε 001.0=ε
Timp de calcul 0.451s 0.600s 0.611s
Pe baza celor prezentate putem trage concluzii legate de comportarea
algoritmilor de tip „quasi-newton”:
Pe fiecare iteraţie este necesară numai valoarea funcţiei gradient, fără a fi necesară evaluarea şi inversarea matricei hessian sau soluţionarea unui sistem de ecuaţii algebrice de mari dimensiuni.
Dacă funcţia obiectiv este o funcţie pătratică atunci convergenţa şi soluţia de optim se va face în cel mult n paşi (n reprezintă dimensiunea vectorului parametrilor).
Indiferent de iniţializarea nx ℜ∈ , metoda converge către un punct de minim (eventual minim local) al criteriului de calitate impus.
Cum matricele de deflexie sunt obţinute prin acumularea informaţiilor obţinute în precedentele iteraţii, metodele quasi-newton sunt mai sensibile la erori de calcul decât metoda Newton. Pentru aceste metode de calculul, vectorul gradient este necesar să fie foarte exact.
Viteza de convergenţă a algoritmului este puternic dependentă de funcţia obiectiv şi de punctul de iniţializare al căutării.
Algoritmul este lent într-o primă fază de lucru. Pe măsură ce ne apropriem de punctul de minim (caz în care aproximarea pătratică este
tot mai bună iar 1−≈ HDk ), algoritmul devine tot mai eficient. Dacă
funcţia criteriu este de clasă C2 şi cu hessianul pozitiv definit în punctul optim, pentru ∞→k , convergenţa este superliniară şi deci algoritmul este mai eficient decât un algoritm de tip gradient. Dacă în plus, în vecinătatea punctului optim, hessianul este lipschitzian, algoritmul are o convergenţă asimptotică cuadratică (comportarea este similară cu algoritmul Newton).
Tabelul 12.4. Timpul de procesare în funcţie de precizie
103
Este posibil, ca din cauza erorilor de calcul, pentru un k, matricea de
deflexie Dk să fie singulară. În cest caz, se reiniţializează Dk=I şi
algoritmul demarează cu o procedură de tip gradient. Pentru a evita
testarea permanentă a condiţiei 0>kD , se poate ca după n iteraţii să
reiniţializăm cu D1 = I.
Metoda poate fi privită ca o tehnică de gradient în care matricea este
modificată pe fiecare iteraţie în ideea sfericizării hipersuprafeţelor de
izonivel. Din acest motiv, tehnica prezentată mai este cunoscută ca
„metodă de metrică variabilă”.
3. Chestiuni de studiat
Se consideră funcţia (forma funcţiei va fi impusă de conducătorul lucrării) a
cărei formă este prezentată în tabelul de mai jos:
Nr. Crt. Funcţia Punct de
iniţializarePunct deminim
Valoare de
minim
1. 2 21 1 1 2 26 2 2 2x x x x x⋅ + ⋅ − ⋅ ⋅ + ⋅ [ ]1, 1− − [ ]2, 1− − -6
2. ( )221 2 1 2
1 109
x x x x+ + ⋅ + − [ ]1, 1− − [ ]5,0.5 7.5
3. ( ) ( )2 21 1 21 100x x x− + ⋅ − [ ]3, 4 [ ]1,1 0
4. ( ) ( )2 21 25 3 5x x⋅ − + − [ ]0,0 [ ]3,5 0
5. 2 21 1 2 2x x x x− ⋅ + [ ]1, 2 [ ]0,0 0
6. 2 21 2 1 29 16 90 128x x x x⋅ + ⋅ − ⋅ − ⋅ [ ]0,3 [ ]5,4 -481
7. 2 21 2 1 2
1 2
2 2 24 6⋅ + ⋅ + ⋅ ⋅ −
− ⋅ − ⋅x x x x
x x [ ]1,1
1 4,3 3⎡ ⎤⎢ ⎥⎣ ⎦
143
−
104
8. 2 21 1 2 2 1 22x x x x x x− ⋅ + − ⋅ + [ ]3,5 [ ]1,0 -1
9. 2
1 1 2 2 1
2
5 4 1612⋅ + ⋅ ⋅ + − ⋅ −
− ⋅x x x x x
x [ ]1,1 [ ]4,14− -152
10. 2 21 2 1 2 1
2
2 2 118⋅ + ⋅ + ⋅ − ⋅ −
− ⋅x x x x x
x [ ]3, 5− − [ ]2,3 -23
11. 2 21 2 1 1 2 22 2x x x x x x− + ⋅ + ⋅ ⋅ + [ ]1,1
31,2
⎡ ⎤−⎢ ⎥⎣ ⎦ -1.25
12. 2 21 2 1 2x x x x+ + ⋅ [ ]1,1 [ ]0,0 0
13. 2 21 216x x+ ⋅ [ ]2,2 [ ]0,0 0
14. ( ) ( )2 21 1 21 x x x− + − [ ]5, 8− − [ ]1,1 0
a) Să se traseze graficul funcţiei considerate împreună cu curbele de izonivel.
b) Să se evalueze eventualul punct de minim utilizând algoritmul propus.
c) Să se reprezinte grafic evoluţia în procesul de căutare.
105
Lucrare de laborator Metoda direcţiilor admisibile de coborâre
(cazul liniar)
1. Scopul lucrării
Scopul principal al acestei lucrări este ca prin implementarea şi
rularea unui algoritm bazat pe metoda direcţiilor admisibile de coborâre să
prezinte probleme de optimizare cu restricţii impuse.
2. Prezentarea lucrării
În cazul problemelor de optim fără restricţii, procedurile clasice de
căutare se bazează pe stabilirea unei direcţii de coborâre şi evoluţia în lungul
acestei direcţii cu un pas de căutare astfel încât să realizăm condiţia de
modificare în sens descrescător a valorilor funcţiei obiectiv. În general,
aplicarea unei astfel de metode nu este aplicabilă în cazul minimelor cu
restricţii deoarece evoluţia în lungul direcţiei de coborâre conduce într-o
zonă în care punctele de lucru nu aparţin domeniului valorilor admise. Prin
urmare ideea de „direcţie de coborâre” îşi pierde aplicabilitatea în cazul
problemelor de minim în prezenţa restricţiilor. Considerăm un punct de
lucru curent Uxk ∈ în care nRU ⊂ reprezintă mulţimea punctelor ce
satisfac restricţiile impuse.
Direcţia de căutare dk conformă cerinţelor problemei necesită să
satisfacă simultan două condiţii:
• pentru un )()(,0 kk xfdxf <⋅+> λλ ;
• pentru λ impus anterior, Udx k ∈⋅+ λ .
13
106
O astfel de direcţie, evident dependentă de xk, poartă numele de
direcţie de coborâre admisibilă.
Pentru început ne punem problema modalităţii construcţiei direcţiilor
admisibile de coborâre.
Considerăm problema minimizării funcţiei RRxf n →:)( , în
restricţia nk RUx ⊂∈ . Pentru un punct Ux∈ , d constituie o direcţie de
coborâre dacă pentru un 0>δ arbitrar de mic, )()( xfdxf <⋅+ λ şi
Udx ∈⋅+ λ pentru ),0( δλ∈∀ .
Pentru început, vom analiza cazul particular al restricţiilor liniare.
Problema generală este de forma:
minimizează f(x)
în restricţiile bxA ≤⋅ (1)
exE =⋅
în care nlmnm RERbRA ××× ∈∈∈ ,, 1 şi 1×∈ lRe . Prin urmare mulţimea de
admisibilitate va fi { , }U x A x b E x e= ⋅ ≤ ⋅ = . Pe motive lesne de înţeles, o
astfel de mulţime este definită într-o serie de lucrări ca „mulţime
poliedrală”).
Legat de modul în care sunt satisfăcute restricţiile de tip inegalitate,
remarcăm că pentru un x fixat este posibil ca relaţia la pasul i să fie
satisfăcută ca egalitate:
ininii bxaxaxa =⋅++⋅+⋅ ...2211 . (2)
Restricţiile de acest gen sunt restricţii active. În mod evident
restricţiile Ex = e, sunt permanent active.
107
În continuare, vom prezenta un rezultat esenţial în procedura de
construcţie a direcţiilor admisibile de căutare. Vom considera problema
generală şi un punct Ux∈ pentru care, printr-o reordonare a matricei A,
restricţiile sunt satisfăcute în forma 2211 , bxAbxA <⋅=⋅ şi exE =⋅ .
Matricea A a fost reordonată şi partiţionată încât
⎟⎟⎠
⎞⎜⎜⎝
⎛=
2
1
AA
A
să fie puse în evidenţă direct restricţiile active (A1, E) pentru punctul
considerat. În aceste condiţii, d reprezintă o direcţie admisibilă de căutare
pentru Sx∈ dacă şi numai dacă 01 ≤⋅ dA şi 0=⋅ dE .
Deoarece Ux∈ avem satisfăcute:
.
,
22
11
exEbxAbxA
=⋅<⋅=⋅
(3)
Pentru λ suficient de mic este necesar ca:
edExEdxE
bdAxAdxAbdAxAdxA
=⋅⋅+⋅=⋅+≤⋅⋅+=⋅+≤⋅⋅+=⋅+
λλλλλλ
)()()(
2222
1111
(4)
Pentru primul set de inegalităţi, cum A1x=b1 rezultă 01 ≤⋅⋅ dAλ şi
deci 01 ≤⋅ dA . Cel de al doilea set de inecuaţii impune
xAbdA ⋅−≤⋅⋅ 222λ şi cum 022 >⋅− xAb putem stabili λ suficient de
mic care să satisfacă condiţia. În sfârşit cu exE =⋅ rezultă necesitatea cu
0=⋅ dE .
Dacă în plus, direcţia admisibilă de căutare satisface 0)( <⋅∇ Τ dxf ,
ea se transformă în direcţie „admisibilă de coborâre”.
108
Vom nuanţa elementele mai sus definite pe următorul exemplu. Se
consideră următoarea problemă de optimizare:
minimizeazǎ (x1 - 6)2 + (x2 - 2)2
în restricţiile 42 21 ≤⋅+− xx
1223 21 ≤⋅+⋅ xx
01 ≤− x
02 ≤− x
Considerăm punctul curent 3,2 21 == xx pentru care:
4322 =⋅+− , 123223 =⋅+⋅ , - 2 < 0 şi -3 < 0.
Rezultǎ matricele A1, A2 şi E în forma:
⎟⎟⎠
⎞⎜⎜⎝
⎛−=
2321
1A , ⎟⎟⎠
⎞⎜⎜⎝
⎛−
−=
1001
2A şi E = 0.
Reprezentarea graficǎ a acestei probleme este urmǎtoarea:
Figura 27. Direcţia admisibilǎ de coborâre
109
Zona haşuratǎ reprezintǎ evoluţia de coborâre admisǎ, iar funcţia f(x)
are o evoluţie admisǎ doar în acea arie.
Considerǎm direcţia ⎟⎟⎠
⎞⎜⎜⎝
⎛=
2
1
dd
d . Condiţia ca sǎ fie direcţie admisibilǎ
este: ⎟⎟⎠
⎞⎜⎜⎝
⎛=≤⎟⎟
⎠
⎞⎜⎜⎝
⎛⋅⎟⎟⎠
⎞⎜⎜⎝
⎛−≤⋅
2
1
2
11 ;0
2321
;0dd
ddd
dA .
Gradientul funcţiei este:
⎟⎟⎠
⎞⎜⎜⎝
⎛−=∇⎟⎟
⎠
⎞⎜⎜⎝
⎛−⋅−⋅
=
⎟⎟⎟⎟⎟
⎠
⎞
⎜⎜⎜⎜⎜
⎝
⎛
∂∂∂∂
=∇=
=
28
)(;)2(2)6(2
)(3
2
2
1
2
1
2
1
x
xxfxx
xfxf
xf .
Pentru ca direcţia consideratǎ sǎ fie direcţie de coborâre, este necesar ca:
028;0)( 21 <⋅+⋅−≤⋅∇ dddxf xT .
Rezultatele anterior prezentate indică condiţiile necesare şi suficiente
ca direcţia d asociată unui punct Ux∈ să fie direcţie admisibilă de
coborâre şi anume
0)(,0,01 <⋅∇=⋅≤⋅ Τ dxfdEdA . (5)
O asemenea condiţionare nu permite stabilirea unei direcţii
admisibile de coborâre unice (a se vedea exemplu prezentat). În ipoteza
existenţei mai multor soluţii se pune problema stabilirii soluţiei care asigură
eficienţa maximă în coborâre (categoric, şi în acest caz ca şi în cazul
tehnicilor de gradient evaluarea este locală), pentru o cât mai rapidă
coborâre este necesar ca dxf ⋅∇ Τ )( să aibă o valoare cât mai mică. Fără
restricţionări suplimentare, dacă există direcţii d cu 0)( * <⋅∇ Τ dxf atunci
110
soluţia problemei va fi ∞− , cu vectorul soluţiei *dµ ⋅ caracterizat de un
∞→µ . Prin urmare problema implică o normare care să limiteze valorile
vectorului d.
Soluţionarea efectivă poate fi realizată cu diverse condiţii de normare:
P1. minimizează dxf )(Τ∇
în restricţiile 01 ≤dA
nidi
Ed
,1,11
0
∈≤≤−
=
P2. minimizează dxf )(Τ∇
în restricţiile 01 ≤dA
1)(
0−≥∇
=Τ dxf
Ed
P3. minimizează dxf )(Τ∇
în restricţiile 01 ≤dA
1
0≤
=Τdd
Ed
Primele două probleme sunt probleme de programare liniară iar ultima
este o problemă de programare pătratică. Considerăm că problema
programării liniare şi a programării pătratice sunt cunoscute atât sub aspect
teoretic cât şi implementativ. Din acest motiv ne vom rezuma la o serie de
observaţii legate strict de contextul general al problemei.
Este evident că d = 0 este o soluţie admisibilă pentru oricare din
problemele P1, P2 sau P3. Pentru o astfel de alegere, valoarea criteriului de
calitate dxf ⋅∇ Τ )( este nulă. Prin urmare, dacă există o soluţie admisibilă
111
care minimizează criteriul atunci valoarea de minim va fi mai mică sau cel
mult egală cu zero. În ipoteza că direcţia d obţinută ca soluţie de minim
asigură 0)( <⋅∇ Τ dxf atunci direcţia obţinută constituie direcţia admisibilă
de căutare şi algoritmul continuă. Dacă soluţia problemei d stabileşte
valoarea indicelui 0)( =⋅∇ Τ dxf atunci x, constituie un punct
Kuhn-Tucker şi algoritmul se consideră încheiat. Vom prezenta în
continuare poate cel mai cunoscut algoritm ce utilizează căutarea pe direcţii
admisibile de coborâre şi anume algoritmul Zoutendijk.
Înainte de a trece la prezentarea algoritmului, sunt necesare câteva
precizări legate de alegerea lungimii pasului de căutare. Am prezentat
anterior condiţii necesare astfel ca d să fie o direcţie admisibilă de coborâre:
numai dacă 01 ≤⋅ dA , 0=⋅ dE şi 0)( <⋅∇ Τ dxf . Singura inegalitate care
condiţionează valorile lui 0>λ este:
2222 )( bdAxAdxA ≤⋅⋅+⋅=⋅+ λλ , (6)
în condiţiile bxA <2 . Prin urmare
xAbdA ⋅−≤⋅⋅ 222λ . (7)
Notăm în continuare:
dAcxAbl ⋅=⋅−= 222 , . (8)
În aceste condiţii, limitarea valorilor R∈λ pentru a asigura condiţia
de coborâre va fi maxλλ ≤ unde maxλ se obţine în forma:
⎪⎩
⎪⎨
⎧
≤⎭⎬⎫
⎩⎨⎧
>
≤∞
=0,0min
0,
max cdacăccl
cdacă
ii
iλ (9)
112
Prezentăm în continuare schiţa algoritmului Zontendijk pentru cazul
restricţiilor de tip egalitate:
Etapa de iniţializare. Trebuie stabilit punctul de iniţializare a căutării
),( 111 dxEbxAUx =⋅≤⋅∈ . Considerăm k=1 şi trecem la etapa de bază.
Etapa de bază. Pasul 1. Considerăm că printr-o reordonare a inecuaţiilor
restricţii obţinem partiţiile A1, A2 astfel ca exEbxA kk =⋅=⋅ ,11 şi
22 bxA k <⋅ . Se rezolvă problema de programare liniară:
minimizează dxf ⋅Τ )(
în restricţiile 01 ≤⋅ dA
njd
dE
j ,1,11
0
∈≤≤−
=⋅
Dacă 0)( =⋅∇ Τkk dxf atunci xk este un punct Kuhn-Tucker şi
algoritmul se opreşte, dacă nu se trece a pasul 2.
Pasul 2. Se determină lungimea pasului de căutare
)(minargmax0
kkk dxf ⋅+=≤<
λλλλ
.
Se determină noul punct de lucru kkkk dxx ⋅+=+ λ1 . Pentru acest punct
determinăm restricţiile active A1, A2. Iterăm 1+→ kk şi reluăm pasul 1.
În continuare vom prezenta aplicarea metodei prezentate pe un
exemplu concret.
Fie următoarea funcţie obiectiv:
2 2( , ) 2 2 2 4 6f x y x y x y x y= ⋅ + ⋅ − ⋅ ⋅ − ⋅ − ⋅
113
Se cere determinarea punctului de minim care asigură următoarele restricţii:
25 5
0
x yx yx
y o
+ ≤⎧⎪ + ⋅ ≤⎪⎨ ≥⎪⎪ ≥⎩
Utilizând subrutina dac.m care implementează algoritmul propus,
obţinem următorul rezultat dintr-o iniţializare 0 [0 ; 0]x = .
x =
0 0.8333 1.1290
0 0.8333 0.7742
-1 -0.5 0 0.5 1-1
-0.5
0
0.5
1
Initializarea cautarii
Punct final
Figura 28. Graficul în procesul de căutare
114
3. Chestiuni de studiat
Se consideră funcţia obiectiv împreună cu restricţiile corespunzătoare una
din funcţiile prezentate în tabelul de mai jos (funcţia va fi impusă de
conducătorul lucrării):
Nr.
crt. Funcţia Restricţii
Pct.
iniţ. f(x0) Minim
Valoare
de
minim
1.
2 21 2
1 2 1
2
2 22 46
⋅ + ⋅ ++ ⋅ ⋅ − ⋅ −− ⋅
x xx x xx
1 2
1
2
2 200
x xxx
− + ⋅ ≤≥≥
0.5,0.5⎡ ⎤⎢ ⎥⎣ ⎦
-3.5
1 ,356
⎡ ⎤⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎣ ⎦
4.16
2. 1 21 2
x xe x x− − − 1 2
1
2
100
x xxx
+ ≤≥≥
[ ]0,0
1 [ ]0,1 -0.6
3.
2 21 2 1
2
42+ − ⋅ −
− ⋅x x x
x
1 2
1 2
1
2
2 42 600
x xx xxx
⋅ + ≤+ ⋅ ≤≥≥
[ ]1,1 -4 1.6,0.8⎡ ⎤⎢ ⎥⎣ ⎦
4.8
4.
2 22 1
1 2
0.5 0.52 5
⋅ + ⋅ −− − ⋅ +
x xx x
1 2
1 2
1
2
2 3 64 500
x xx xxx
⋅ + ⋅ ≤+ ⋅ ≤≥≥
[ ]0,0
5
13 ,171817
⎡ ⎤⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎣ ⎦
10134
5.
21 1
22 2
2 0.2
3 0.2
− ⋅ + ⋅ −
− ⋅ + ⋅
x x
x x
1 2
1 2
1
2
2 3 132 10
00
x xx x
xx
⋅ + ⋅ ≤⋅ + ≤≥≥
[ ]1,1 -4.6 [ ]2,3 -10.4
6.
2 21 2 1
2
2.45.6+ − ⋅ −
− ⋅x x x
x
1 2
1 2
1 2
1
2
2 3 3 03 0
2 4 000
x xx x
x xxx
− ⋅ − ⋅ − ≤+ − ≤⋅ − − ≤≥≥
[ ]0,1 -4.6
1.2,1.8⎡ ⎤⎢ ⎥⎣ ⎦
-7.88
115
a) Să se traseze graficul funcţiei considerate împreună cu curbele de izonivel.
b) Să se evalueze eventualul punct de minim utilizând algoritmul propus.
c) Să se reprezinte grafic evoluţia în procesul de căutare.
Lucrare de laborator Metoda gradientului proiectat
(cazul liniar)
1. Scopul lucrării Scopul principal al acestei lucrări este ca prin implementarea şi rularea unui
algoritm bazat pe metoda gradienţilor proiectaţi, pentru cazul liniar.
2. Prezentarea lucrării Considerăm problema generală de optimizare în prezenta restricţiilor:
minimizează f (x)
în restricţia nRUx ⊂∈ , (1)
unde U reprezintă mulţimea punctelor admisibile, iar f (x) funcţia continuă
cu derivate parţiale continue pe U ( )()( 1 UCxf ∈ ). În cadrul metodei celei
mai rapide coborâri, căutarea minimului se face prin construcţia iterativă
)(1 kkkk xfxx Τ+ ∇⋅−= λ . (2)
În cazul în care direcţia de antigradient nu este o direcţie admisibilă în raport cu mulţimea U, se propune ca direcţie de căutare proiecţia lui xk + 1 pe
mulţimea U.
,...2,1))((1 =∇⋅−= Τ+ kxfxPx kkkuk λ (3)
14
116
în care 0>kλ . Dacă mulţimea U este convexă, succesiunea de puncte de
căutare {xk} este perfect definită. În cazul particular, U = Rn metoda
coincide cu metoda celei mai rapide coborâri.
Dacă mulţimea U este convexă şi dacă )()( 1 UCxf ∈ notând U*
mulţimea punctelor de minim pentru funcţia f(x) cu Ux∈ atunci
))(( *** xfxPx uΤ∇⋅−= λ (4)
În baza acestui rezultat şi în ipoteza că în cadrul procesului iterativ
(2), xk+1=xk atunci procesul iterativ se opreşte şi punctul xk satisface
condiţiile necesare de optim. Dacă în plus, f(x) este o funcţie convexă,
atunci xk = x* reprezintă soluţia de optim.
În funcţie de modul de alegere a lungimii pasului de căutare kλ ,
distingem mai multe proceduri de optimizare bazate pe gradienţi proiectaţi.
Prezentăm în continuare cele mai utilizate proceduri de alegere a pasului de
căutare:
1) Definim funcţia scalară pe iteraţia k:
,...2,1))((()( =∇⋅−= Τ kxfxPIf kkuk λλ (5)
unde 0>λ . Lungimea pasului de căutare se determină soluţionând
problema de minim unidimensional
0)(inf)( *
0>==
≥ kkkkk fff λλλλ
. (6)
Evident, pentru cazul U = Rn procedura iterativă (5) cu lungimea
pasului de căutare kλ determinat de (6), coincide cu metoda celei mai rapide
coborâri.
117
2) Există posibilitatea determinării prin testare succesivă a unei
valori kλ care asigură condiţia de monotonie, )()( 1 kk xfxf <+ .
3) Dacă funcţia f(x) continuă, cu derivate parţiale continue, asigură
condiţia Lipschitz pentru funcţia gradient cu constanta L, putem alege
pentru lungimea pasului de căutare valoarea
)2(
20 0 Σ+≤≤Σ<
Lkλ , (7)
unde oΣΣ, sunt parametri de metodă.
4) Este posibil ca lungimea pasului de căutare kλ să fie predefinit cu
asigurarea următoarelor restricţii:
∑∞
=
∞==>1
,...,2,10K
kk k λλ şi ∑∞
=
∞<1
2
Kkλ , (8)
ca de exemplu 1
1+
=kkλ .
Pentru accelerarea procesului de căutare, în locul iteraţiei (3) se
poate utiliza
0,10)1())((
)))(((1
>≤<⋅−+∇−=
=−∇−+=Τ
Τ+
kkkkkkkuk
kkkkukkk
xxfxP
xxfxPxx
λββλβ
λβ (9)
în care kλ şi kβ pot fi aleşi în diverse moduri.
Reamintim că indiferent de modul de alegere a lungimii pasului de
căutare, pe fiecare iteraţie este necesar să proiectăm punctul rezultat pe
mulţimea restricţiilor. Altfel spus, în cadrul fiecărei iteraţii trebuie
soluţionată o problemă de forma:
2
min ( ( )) ,k k kx x f x x Uλ Τ− − ⋅∇ ∀ ∈ (10)
118
Din acest motiv, metoda poate fi aplicată numai pentru cazurile în care restricţiile impuse determină o mulţime U pentru care proiecţia se obţine relativ simplu.
În continuare, vom prezenta un caz, pentru care proiecţia gradientului pe mulţimea generată de restricţii se obţine relativ simplu (metoda gradientului proiectat).
Considerăm cazul minimizării unei funcţii nRxf :)( , pentru care
restricţiile de tip egalitate sau inegalitate sunt liniare:
minimizează f(x) în restricţiile bxA ≤⋅ (11)
exE =⋅
unde nmRA ×∈ , 1×∈ mRb , nlRE ×∈ şi 1×∈ lRe . Dacă considerăm un punct de
lucru curent nRx∈ şi dacă acest punct este interior mulţimii U, oricare direcţie este o direcţie admisibilă şi o evoluţie pe direcţia de antigradient este posibilă dacă alegem un pas suficient de mic pentru a nu părăsi frontiera mulţimii U. În cazul în care punctul aparţine frontierei mulţimii U, evoluţia pe sensul impus de antigradient impune în general evoluţia în puncte ce nu satisfac restricţiile impuse. Ideea metodei constă în proiectarea direcţiei de antigradient pe frontiera domeniului şi organizarea căutării pe o astfel de direcţie.
Fie punctul Ux∈ pentru care putem stabili partiţii corespunzătoare astfel încât A1x = b, A2x <b şi Ex = e. Evident, partiţiile A1, A2, b1 şi b2 se obţin printr-o reordonare corespunzătoare a liniilor matricelor A şi b.
Restricţiile active pentru punctul nRx∈ considerat, fixează o matrice
⎟⎟⎠
⎞⎜⎜⎝
⎛=
EA
M 1
119
pe care o presupunem de rang complet. Introducem matricea P =I -
MT(MMT)-1M. Evident, matricea P este simetrică şi involutivă şi prin urmare
P este o matrice proiector. Aplicând operatorul P oricărui vector, obţinem
proiecţia acestui vector în subspaţiul restricţiilor active. Pe ideea generală,
vom propune ca direcţie de căutare, proiecţia antigradientului în subspaţiul
restricţiilor active, deci )(xfPd ∇⋅−= .
Vom arăta că în cazul în care 0)( ≠∇⋅ xfP , direcţia )(xfPd ∇⋅−=
este o direcţie admisibilă de căutare. În adevăr:
0)()()()()()( 2 <∇⋅−=∇⋅⋅−∇=∇⋅⋅−∇=⋅∇ ΤΤΤΤ xfPxfPPxfxfPxfdxf , (12)
şi deci d este o direcţie de coborâre. Pe de altă parte,
0)( =∇⋅⋅−=⋅ xfPMdM şi deci A1d = 0 , Ed = 0 care reprezintă condiţia
ca d să fie o direcţie admisă.
În cazul în care 0)( =∇ xfP , obţinem:
)()()()())(()( 11 xfMMMMxfxfMMMMIxfP ∇⋅⋅⋅⋅−∇=∇⋅⋅⋅⋅−=∇ −ΤΤ−ΤΤ
Vom nota în continuare
⎟⎟⎠
⎞⎜⎜⎝
⎛==∇⋅⋅− −Τ
vu
wxfMMM )())( 1 .
În aceste condiţii,
0)()()( 11 =⋅+⋅−∇=⎟⎟
⎠
⎞⎜⎜⎝
⎛⋅⎟⎟
⎠
⎞⎜⎜⎝
⎛+∇=∇⋅ ΤΤ
Τ
vEuAxfvu
EA
xfxfP . (13)
Facem observaţia că A1 şi E corespund gradienţilor restricţiilor active
şi prin urmare dacă toate componentele lui u asigură 0≥u atunci punctul x
este un punct Kuhn-Tucker.
120
Dacă vectorul u nu satisface condiţia 0≥u , există cel puţin o
componentă 0<ju . Considerăm că în matricea A1 eliminăm linia j obţinând
matricea Â1 şi notăm ⎟⎟⎠
⎞⎜⎜⎝
⎛=
EA
P 1ˆ . Vom arăta că 0)(ˆ ≠∇⋅ xfP .
Presupunem prin absurd că 0)(ˆ =∇⋅ xfP şi deci
0ˆ)(
)(ˆ)ˆˆ(ˆ)(
)()ˆ)ˆˆ(ˆ(1
1
=⋅+∇=
=∇⋅⋅⋅⋅−∇=
=∇⋅⋅⋅⋅−
Τ
−ΤΤ
−ΤΤ
WMxf
xfMMMMxf
xfMMMMI
(13)
în care WxfMMM =∇⋅⋅− −Τ )(ˆ)ˆˆ( 1 .
Pe de altă parte, jjurWMvEuA ˆˆˆ1 +=+ ΤΤΤ , unde rj reprezintă un
vector corespunzător liniei j din A1. Prin urmare,
jj urWMxf ΤΤ ++∇= ˆˆ)(0 . (14)
Scăzând egalitatea (13) şi (14), obţinem:
0)ˆ(ˆ =+− ΤΤjj urWWM . (15)
Cu 0≠in , egalitatea contrazice ipoteza că M este de rang complet.
Prin urmare, 0)(ˆ ≠∇⋅ xfP . Prin modul de definire, matricea P este şi ea o
matrice de proiecţie şi considerând )(ˆ xfPd ∇⋅−= , vom arăta că d
constituie o direcţie admisibilă de coborâre.
Deoarece P̂ este proiector, )(ˆ xfPd ∇⋅−= este o direcţie de coborâre.
Pentru a arăta că d este o direcţie admisibilă, este suficient să arătăm că
01 ≤⋅ dA şi Ed = 0. Ţinând cont că 0ˆˆ =⋅ PM , rezultă:
0)(ˆˆˆ1 =∇⋅−=⋅⎟⎟⎠
⎞⎜⎜⎝
⎛xfPMd
EA , (16)
121
şi deci 0ˆ1 =⋅ dA , Ed = 0. Rămâne să arătăm că 0≤⋅Τ drj .
Înmulţind egalitatea (4.66) cu Prjˆ⋅ , obţinem:
ΤΤΤ ⋅⋅⋅+⋅−=⋅+⋅⋅⋅+∇⋅⋅= jjjjjjjj rPrudrurWMPrxfPr ˆ)ˆˆ(ˆ)(ˆ0 . (17)
Relaţia a fost obţinută ţinând cont că 0ˆˆ =⋅ ΤMP .
nl
nm
RE
RA×
×
∈
∈1 , nlmRMEA ×+∈=⎟⎟⎠
⎞⎜⎜⎝
⎛ )(1 , (18)
⎥⎦
⎤⎢⎣
⎡=∈∇⋅⋅−= ×+
×
−Τ
+×+ vu
RxfMMMw nlm
xnlmln
)(1
)()()()( , (19)
[ ] [ ] [ ] jj
l
m
j
l
m
j
ur
v
vu
u
u
EA
v
vu
u
u
EA ΤΤΤΤΤ +
⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥
⎦
⎤
⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢
⎣
⎡
⋅=
⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥
⎦
⎤
⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢
⎣
⎡
⋅
.
.
.
.
.
.
.
.
.
ˆˆ
.
.
.
.
.
.
.
.
.
1
1
1
1
1
1. (20)
Cum uj < 0 şi 0ˆ >⋅⋅ Τjj rPr ( P̂ proiector şi deci matrice pozitiv definită),
rezultă că 0≤⋅ drj .
Pe baza celor stabilite, prezentăm un algoritm bazat pe metoda gradientului proiectat pe restricţii liniare. În continuare, prezentăm un algoritm (Rosen) pentru minimizarea funcţiilor scalare de mai multe
variabile RRxf n →:)( , în cazul restricţiilor poliedrale exEbxA =⋅≤⋅ , .
122
Etapa de iniţializare. Se impune nRx ∈1 pentru a satisface restricţiile
impuse bxA ≤⋅ 1 şi exE =⋅ 1 . Reordonăm restricţiile astfel încât putem
fixa partiţiile ⎟⎟⎠
⎞⎜⎜⎝
⎛=⎟⎟
⎠
⎞⎜⎜⎝
⎛=
2
1
2
1 ,bb
bAA
A încât 111 bxA =⋅ şi 212 bxA <⋅ . Facem k =
1 şi trecem la etapa de bază.
Etapa de bază. Pasul 1. Fie ⎟⎟⎠
⎞⎜⎜⎝
⎛=
EA
M 1 . Dacă M nu are nici un element,
atunci P = In, în caz contrar, MMMMP ⋅⋅⋅−= −ΤΤ 1)(1 . Calculăm
)( kxfPd ∇⋅−= . Dacă 0≠kd , trecem la pasul 2.
Dacă 0=kd , calculăm )()( 1kxfMMMW ∇⋅⋅−= −Τ . Fie ⎟⎟
⎠
⎞⎜⎜⎝
⎛=
vu
W . Dacă
0≥u ne oprim, iar xk este un punct Kuhn-Tucker. Dacă nu, 0≥u atunci
alege, 0<ju , eliminăm din A1 linia j-a şi reluăm pasul 1.
Pasul 2. Se determină lungimea optimă a pasului de căutare kλ soluţionând
problema:
minimalizează )( kk dxf ⋅+ λ
max0 λλ ≤≤
unde maxλ se determină ca şi în cazul metodei direcţiilor admisibile de
căutare. Construim kkkk dxx ⋅+=+ λ1 , stabilim partiţiile A1, A2 ale lui A şi
b1, b2 ale lui b încât 111 bxA k =⋅ + , 112 bxA k <⋅ + .
Facem 1+→ kk şi reluăm pasul 1.
3. Chestiuni de studiat
Se consideră funcţia obiectiv împreună cu restricţiile corespunzătoare una
din funcţiile prezentate în tabelul următor (funcţia va fi impusă de
conducătorul lucrării):
123
Nr.
crt. Funcţia Restricţii
Pct.
iniţ. f(x0) Minim
Valoare
de
minim
1.
2 21 2
1 2 1
2
2 22 46
⋅ + ⋅ ++ ⋅ ⋅ − ⋅ −− ⋅
x xx x xx
1 2
1
2
2 200
x xxx
− + ⋅ ≤≥≥
0.5,0.5⎡ ⎤⎢ ⎥⎣ ⎦
-3.5
1 ,356
⎡ ⎤⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎣ ⎦
4.16
2. 1 21 2
x xe x x− − − 1 2
1
2
100
x xxx
+ ≤≥≥
[ ]0,0
1 [ ]0,1 -0.6
3.
2 21 2 1
2
42+ − ⋅ −
− ⋅x x x
x
1 2
1 2
1
2
2 42 600
x xx xxx
⋅ + ≤+ ⋅ ≤≥≥
[ ]1,1 -4 1.6,0.8⎡ ⎤⎢ ⎥⎣ ⎦
4.8
4.
2 22 1
1 2
0.5 0.52 5
⋅ + ⋅ −− − ⋅ +
x xx x
1 2
1 2
1
2
2 3 64 500
x xx xxx
⋅ + ⋅ ≤+ ⋅ ≤≥≥
[ ]0,0
5
13 ,171817
⎡ ⎤⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎣ ⎦
10134
5.
21 1
22 2
2 0.2
3 0.2
− ⋅ + ⋅ −
− ⋅ + ⋅
x x
x x
1 2
1 2
1
2
2 3 132 10
00
x xx x
xx
⋅ + ⋅ ≤⋅ + ≤≥≥
[ ]1,1 -4.6 [ ]2,3 -10.4
6.
2 21 2 1
2
2.45.6+ − ⋅ −
− ⋅x x x
x
1 2
1 2
1 2
1
2
2 3 3 03 0
2 4 000
x xx x
x xxx
− ⋅ − ⋅ − ≤+ − ≤⋅ − − ≤≥≥
[ ]0,1 -4.6
1.2,1.8⎡ ⎤⎢ ⎥⎣ ⎦
-7.88
a) Să se traseze graficul funcţiei considerate împreună cu curbele de izonivel.
b) Să se evalueze eventualul punct de minim utilizând algoritmul propus.
c) Să se reprezinte grafic evoluţia în procesul de căutare.