Optimizare fara restrictii

3

OPTIMIZARE MULTIDIMENSIONAL FR RESTRICII

3.1 Introducere n multe situaii practice este necesar determinarea minimului unei funcii obiectiv F(X), fr a avea alte condiii suplimentare impuse variabilelor de decizie x1, x2, , xn, coninute n vectorul X. Determinarea soluiei X* a problemei de optimizare fr restricii impune anularea gradientului funciei obiectiv, adic: ( ) ;0* = XF unde gradientul funciei obiectiv este vectorul derivatelor pariale de ordinul nti, definit de relaia (1.8). Din punct de vedere practic, este puin probabil s rezolvm aceast problem direct, deoarece componentele gradientului pot avea expresii puternic neliniare. n astfel de cazuri putem apela la o procedur numeric de minimizare. De cele mai multe ori procedurile numerice necesit alegerea unui punct de start X0. Apoi, n mod iterativ, printr-o metod oarecare, se caut un punct nou X1, n care funcia obiectiv are valoarea mai mic dect n X0 s.a.m.d. Una dintre cele mai simple metode de minimizare a unei funcii dependente de n variabile, fr restricii, const n selecia aleatoare a unui numr mare de valori pentru vectorul X al variabilelor de decizie i evaluarea funciei obiectiv n fiecare dintre aceste puncte. Punctul X corespunztor valorii minime a lui F(X) este considerat punct de optim, sau de minim, notat X*. Firete c dac dorim o soluie foarte precis, atunci va trebui s avem n vedere un numr foarte mare de valori ale vectorului X. Metodele de acest tip, n care pentru determinarea punctului de optim apelm doar la calculul valorilor funciei obiectiv n anumite puncte, poart numele de metode de ordinul zero. Metode puin mai dificile, dar i mai eficiente, sunt cele care pentru determinarea punctului de optim utilizeaz derivata de ordinul nti a funciei obiectiv pentru determinarea direciei de minimizare. tim c gradientul funciei obiectiv este vectorul a crui direcie este ndreptat n sensul creterii maxime a valorii funciei obiectiv. Din acest motiv cutarea punctului de minim se face n sens opus celui indicat de vectorul gradient. Se pleac tot dintr-un punct de start X0 i se caut un nou punct X1, n care funcia obiectiv are valoare mai mic dect n punctul precedent. Relaia de tip iterativ pentru determinarea lui X1 este: ( );001 XFXX += ncercnd mai multe valori pentru scalarul , vom determina acea valoare notat

(3.1)

(3.2)

Optimizare numeric. Algoritmi i programe n C

56

*, pentru care F(X) are valoarea minim n direcia S0 = -F(X0). Folosind aceast valoare *, cu ajutorul relaiei (3.2) vom determina punctul X1, ce reprezint minimul lui F(X) pe direcia S0 = -F(X0). Deoarece funcia obiectiv poate fi neliniar, n noul punct X1 va trebui s calculm o nou direcie de minimizare S1 = -F(X1) i pe aceast direcie s determinm punctul X2 conform relaiei (3.2) s.a.m.d. Metodele de acest tip, ce folosesc gradientul funciei la determinarea direciei de minimizare, poart numele de metode de ordinul nti. S considerm acum dezvoltarea n serie Taylor pn la termenii de ordinul doi inclusiv, a funciei obiectiv F(X) n punctul de start X0:

( ) ( ) ( ) ( ) ;21 000 XXHXXXFXFXF T ++

unde 0XXX = , iar matricea hessian este definit conform relaiei (1.9). Ecuaia (3.3) este o aproximare a funciei F(X) i va fi cu att mai bun cu ct punctul X0 va fi mai apropiat de X*. Dac derivm ecuaia (3.3) n raport cu X obinem:

( ) ( ) ( ) ;00 XXHXFXF += Deoarece n punctul de optim, gradientul funciei este nul, avem: ( )[ ] ( ) ;010 XFXHX = i innd cont c 0XXX = , rezult: ( )[ ] ( ) ;0100 XFXHXX = n noul punct X, calculat cu ajutorul relaiei (3.6), repetm dezvoltarea n serie Taylor conform (3.3), dup care calculm noul punct X s.a.m.d. Pentru c n vederea determinrii punctului de minim a funciei obiectiv am folosit i derivata de ordinul al doilea a funciei F(X), aceste metode poart numele de metode de ordinul doi. Dei principial, sunt metode mai eficiente dect cele de ordinul nti, dificultatea aplicrii lor const n calculul matricei hessiene. Principala utilitate a acestor metode const n completarea sau dezvoltarea unor metode de optimizare de ordinul nti. 3.2 Metoda general de optimizare n general, metodele numerice de optimizare de tip iterativ, folosesc o relaie de recuren de forma:

;*1 qqqq SXX +=+

n care: q este un indice iterativ;

(3.3)

(3.4)

(3.5)

(3.6)

(3.7)

Optimizare multidimensional fr restricii

57

X vectorul variabilelor de decizie; Sq vectorul direciei de minimizare; *q - o mrime scalar, ce caracterizeaz mrimea saltului din punctul

Xq, n punctul Xq+1, n direcia lui Sq. Utilizarea relaiei (3.7) ntr-un algoritm de optimizare presupune ndeplinirea a

dou cerine. Prima cerin presupune determinarea direciei de cutare Sq, iar a doua, determinarea factorului de salt *q . n concluzie, putem spune c un algoritm general de optimizare a unei funcii de n variabile, fr restricii, conine trei componente majore:

a) determinarea direciei de cutare; b) calculul punctului de optim, de-a lungul direciei de cutare

determinate anterior; c) determinarea momentului n care procedura de optimizare

converge ctre o soluie acceptabil. 3.3 Metode de ordinul zero Metodele de optimizare ce folosesc pentru determinarea soluiei optime numai valorile funciei n anumite puncte, sunt mai uor de programat. Adesea aceste metode pot lucra efectiv, cu funcii discontinue i uneori chiar cu valori discrete ale variabilelor de decizie. Dezavantajul const n numrul mare de evaluri ale funciei obiectiv, pn la determinarea soluiei cu o anumit precizie, impus iniial. Aceste metode sunt recomandate n probleme cu puine variabile de decizie, la care evaluarea funciei obiectiv nu necesit un timp mare de compilare. 3.3.1 Metoda cutrii aleatoare Metoda cutrii aleatoare este considerat una dintre metodele cele mai simplu de implementat pe calculator. Ea const din alegerea n mod aleator a unor vectori X din spaiul de cutare definit de variabilele de decizie. n vederea evitrii unei cutri pe un spaiu mai mare dect este necesar, se pot impune nite limite rezonabile pentru variabilele de decizie, de forma:

.,...,2,1; nibxa iii = Majoritatea compilatoarelor dispun de funcii generatoare de numere aleatoare uniform distribuite ntr-un anumit interval [0, RAND_MAX], unde RAND_MAX este valoarea maxim a numrului aleator astfel generat. Pentru compilatorul de Microsoft Visual C++, versiunea 6.0, valoarea lui RAND_MAX est 32767.

(3.8)


58

Deoarece avem nevoie de numere aleatoare uniform distribuite n intervalul [ai,bi], vom mpri fiecare numr aleator, generat n intervalul [0, RAND_MAX], la RAND_MAX, transformnd irul de numere aleatoare ntregi, uniform distribuite n intervalul [0, RAND_MAX], ntr-un ir de numere aleatoare uniform distribuite n intervalul [0,1]. Dm n continuare programul surs al metodei cutrii aleatoare, scris pentru o funcie dependent de dou variabile, dar adaptarea acestuia pentru cazul unei funcii dependente de mai multe variabile este foarte simpl. Pentru nceput, se va actualiza valoarea variabilei N_VAR din prima linie de program. Apoi, n ciclul for cu N pai, se vor aduga noi linii de comand pentru generarea valorilor aleatoare corespunztoare variabilelor suplimentare. De asemenea trebuie modificat apelul funciei pentru numrul actualizat de variabile de decizie.

#define N_VAR 2 // numarul variabilelor de decizie #define N 1e9 // numarul punctelor generate aleator #define epsilon 1e-6 // factor de convergenta #define N_SUB 1000 // numar de subintervale pe fiecare din // axe // functia obiectiv #define F(x,y) 2*x+y+x*x-y*y+(x-y*y)*(x-y*y) #include #include #include #include void main(void) { int nGrid = N_SUB + 1;// nGrid este numarul punctelor in care // calculam valoarea functiei double a[N_VAR]; // limita inferioara a domeniului de cautare double b[N_VAR]; // limita superioara a domeniului de cautare double x1, x2; // variabile de decizie double x1Min, x2Min; // coordonatele punctului de minim double f, fMin; // valoarea curenta respectiv minima a // functiei obiectiv double delta[N_VAR]; // valoarea incrementului pe axa i FILE *fp1; // pointer la fisierul de rezultate "date.txt" // initializam limitele domeniului de cautare for(int i=0;i


59

/* * Initializam generatorul de numere aleatoare astfel incat * secventa de numere sa difere la fiecare alta rulare. */ srand( (unsigned)time( NULL ) ); // calculam incrementul pe fiecare axa for( i=0; i


60

Dac testm eficiena programului pe funcia (1.1) dependent de dou variabile x1 i x2 i alegem un factor de convergen = 10-6, respectiv un numr de puncte generate aleator N = 109, vom obine rezultatele din tabelul 3.1.

Tabelul 3.1 iteraia x1 x2 fmin

1 1.247291 -7.899716 3674.073475 2 0.600909 0.191351 2.036075

99 0.242012 -1.056856 -0.865700 130 -0.417798 -0.801721 -0.980742

2.857 -0.000916 -0.830409 -1.045035 2.948 -0.328684 -0.658284 -1.060279 8.591 -0.299997 -0.724815 -1.078959

58.749 -0.195624 -0.721763 -1.082217 66.224 -0.094913 -0.848720 -1.085249

198.874 -0.211493 -0.763268 -1.093554 900.212 -0.193182 -0.793786 -1.095140

6.677.346 -0.188299 -0.796228 -1.095207 10.720.203 -0.186468 -0.793786 -1.095271 76.013.808 -0.185858 -0.793786 -1.095274

266.610.230 -0.184637 -0.794397 -1.095274 419.283.639 -0.185858 -0.793176 -1.095275 770.279.485 -0.185274 -0.793176 -1.095275 810.486.829 -0.184637 -0.793786 -1.095275 811.397.574 -0.185247 -0.793786 -1.095275

Prin urmare, dup generarea aleatoare a celor N puncte i calculul valorilor funciei obiectiv n fiecare dintre ele, soluia gsit este:

.095275.1,793786.0185247.0

min* =

= fcuX

3.3.2 Metoda drumului aleator

Cu toate c este mai eficient dect metoda cutrii aleatoare, aceast metod este tot de tip aleator. Fa de metoda precedent, noutatea const n determinarea unui vector S a direciei de minimizare, cu ajutorul numerelor aleatoare. Este necesar specificarea iniial a unui punct de start X0. Paii efectuai n vederea determinrii punctului de minim se fac pe baza relaiei de tip iterativ:

;1 qqq SXX +=+ unde este o mrime scalar, numit factor de salt, iar S este vectorul direcie, ale

(3.9)


61

crui componente se determin n mod aleator. Pentru aceasta, se genereaz aleator componentele vectorului S i se ncearc, pentru o valoare mic a factorului de salt, dac direcia S astfel generat, conduce spre descreterea valorii funciei obiectiv. Dac direcia aleatoare S conduce spre descreterea valorii funciei obiectiv, atunci se mrete valoarea factorului de salt cu ajutorul unui increment, dup care se recalculeaz valoarea funciei obiectiv. Procedeul continu n acest mod, mrind valoarea factorului de salt, pn cnd valoarea actual a funciei devine mai mare dect valoarea precedent. n acest moment se poate considera cu aproximaie, c am ajuns la minimul funciei, pe direcia S. Acum se va genera o nou direcie aleatoare i tot procesul continu n acelai mod. Determinarea fiecrei noi direcii S, de minimizare se va efectua prin ncercri repetate, pn cnd se va gsi un vector S, ale crui componente s aib astfel de valori, nct nlocuite n relaia (3.9), s conduc spre descreterea valorilor funciei obiectiv. Ne impunem iniial i un factor de convergen , iar procedura de cutare se consider ncheiat atunci cnd diferena absolut a valorilor funciei ntre dou iteraii succesive devine mai mic dect , de n ori la rnd. Pentru aceasta ne putem fixa un numr de nMax = 5 iteraii succesive, n care criteriul de convergen este ndeplinit. Prin aceast condiie, evitm oprirea programului atunci cnd procedura de cutare converge foarte lent, sau se ajunge ntr-un optim local. Cu toate aceste precauii, aceast metod determin punctul de optim local i este bine s repornim programul din mai multe puncte de start diferite, pentru a ne asigura c punctul de optim determinat este chiar cel cutat. De asemenea un punct sensibil al procedurii este legat de alegerea valorii factorului de salt i respectiv a valorii incrementului cu care mrim valoarea acestuia.

Tabelul 3.2 Punctul de start Factor de salt Increment Factor de convergen

01x

02x

3.0 3.5 0.001 0.001 1e-5

Tabelul 3.3 Numr de iteraii Punctul de minim Valoarea optim

i min1x min2x

minf 992 -0.191765 -0.796921 -1.095275

n tabelul 3.2 sunt date valorile alese pentru punctul de start, factorul de salt, incrementul, respectiv factorul de convergen. Valorile aleatoare pentru determinarea vectorului direcie de minimizare s-au generat n intervalul [-10,10]. Tabelul 3.3 prezint soluia optim determinat n cazul funciei (1.1), dup efectuarea a 992 iteraii. n cele ce urmeaz, este prezentat programul surs de


62

optimizare, prin metoda drumului aleator, scris n limbaj C. Programul este fcut s determine minimul unei funcii dependente de dou variabile, dar orice utilizator cu puine cunotine de programare poate s-l adapteze uor pentru n variabile.

#include #include #include #include // macrodefinitia functiei obiectiv #define f(x1,x2) 2*x1+x2+(x1*x1-x2*x2)+(x1-x2*x2)*(x1-x2*x2) void main(void) { int i=1; // indice iterativ int n=0; // index int nMax = 5; // coordonatele punctului de start X0 double x1zero = 3.0; double x2zero = 3.5; // coordonatele punctului X1 double x1unu, x2unu; // factorul de convergenta double epsilon = 1e-5; // valoarea initiala a factorului de salt double alfa0 = 0.001; double alfa; // incrementul factorului de salt double increment; // coordonatele vectorului directie double s1, s2; // intervalul de valori pentru S double a = -10.0; double b = 10.0; // valori ale functiei obiectiv double F0, F1, Fprecedent, Factual; double ff[2]; // diferenta absoluta a valorilor functiei double deltaF = 1e2; FILE *fp; // calculam valoarea functiei in X0 F0 = f(x1zero,x2zero);


63

Fprecedent = F0; ff[0] = F0; /* Initializam generatorul de numere aleatoare astfel incat * secventa de numere sa difere la fiecare alta rulare. */ srand( (unsigned)time( NULL ) ); // deschidem fisierul "date_out.txt" fp = fopen( "date_out.txt", "w" ); // generam o directie S aleatoare s1 = rand(); s1 = s1/RAND_MAX; s2 = rand(); s2 = s2/RAND_MAX; alfa = alfa0; increment = alfa0; // calculam coordonatele lui X1 x1unu = x1zero + alfa*s1; x2unu = x2zero + alfa*s2; // calculam valoarea functiei in X1 F1 = f(x1unu,x2unu); Factual = F1; for(;;) { while(Factual>=Fprecedent) { // generam o directie S aleatoare s1 = rand(); s1 = s1/RAND_MAX; s1 = a + s1*(b-a); s2 = rand(); s2 = s2/RAND_MAX; s2 = a + s2*(b-a); // calculam coordonatele lui X1 x1unu = x1zero + alfa*s1; x2unu = x2zero + alfa*s2; // calculam valoarea functiei in X1 F1 = f(x1unu,x2unu); Factual = F1; } printf("iteratia=%d fmin=%lf\n", i, F1); fprintf(fp, "iteratia=%d fmin=%lf\n", i, F1); i++; ff[1] = F1; while(Factual


64

Fprecedent = Factual; alfa = alfa + increment; // calculam coordonatele lui X1 x1unu = x1zero + alfa*s1; x2unu = x2zero + alfa*s2; // calculam valoarea functiei in X1 F1 = f(x1unu,x2unu); } deltaF = fabs(F1-Fprecedent); F0 = Factual; x1zero = x1unu; x2zero = x2unu; alfa = alfa0; if(fabs(ff[0]-ff[1])


65

care se va considera drept variabil 03x i se vor fixa constante celelalte componente .a.m.d. Dup efectuarea celor n pai necesari determinrii minimului de-a lungul fiecrei axe, se consider ncheiat o iteraie. Diferena valorilor absolute ale funciei corespunztoare la dou iteraii consecutive, poate fi folosit drept criteriu de oprire a procedurii de cutare.

// macrodefinitia functiei obiectiv #define f(u,v) 2*u+v+(u*u-v*v)+(u-v*v)*(u-v*v) #include #include void main( void ) { int j = 1; // punctul de start double x1zero, x2zero; // limitele intervalului de cautare

f2 0 2

2

0

2

7

6

6

6

5

5

5

4

4

4

4

33

3

33

3

2

2

2

2

2

2

1

1

1 1

1

1

0

0

0

00

0

1 1

X0

X1

X2

Fig. 3.1 Interpretarea geometric a metodei relaxrii.


66

double a1 = -5.0, b1 = 5.0; double a2 = -4.0, b2 = 4.0; double ls1, ls2, ld1, ld2; double x1, x2; double x1Min, x2Min, fMin,f; // numarul de subintervale la cautarea locala double n = 250.0; double increment; double delta=1.0e3; double epsilon = 1.0e-9; double f_actual, f_precedent; FILE *fp; // deschidem fisierul "date_out.txt" fp = fopen( "date_out.txt", "w" ); printf("\n\n METODA RELAXARII (MULTI-GRID)\n"); fprintf(fp, "\n\n METODA RELAXARII (MULTI-GRID)\n"); printf(" ========================================\n\n"); fprintf(fp, " ========================================\n\n"); ls1 = a1; ld1 = b1; ls2 = a2; ld2 = b2; printf(" Limitele intervalului de cautare:\n"); fprintf(fp, " Limitele intervalului de cautare:\n"); printf("\n ls1=%lf ld1=%lf\n ls2=%lf ld2=%lf\n\n", ls1, ld1, ls2, ld2); fprintf(fp, "\n ls1=%lf ld1=%lf\n ls2=%lf ld2=%lf\n\n", ls1, ld1, ls2, ld2); // calculul functieiin punctul de start x1zero = 3.0; x2zero = 3.5; printf(" Coordonatele punctului de start:\n\n"); fprintf(fp, " Coordonatele punctului de start:\n\n"); printf(" x1zero=%lf x2zero=%lf\n\n", x1zero, x2zero); fprintf(fp, " x1zero=%lf x2zero=%lf\n\n", x1zero, x2zero); f_precedent = f(x1zero,x2zero); printf(" F(X0)=%lf\n\n", f_precedent); fprintf(fp, " F(X0)=%lf\n\n", f_precedent); fMin = f_precedent; // incepem din punctul de start


67

x2 = x2zero; do { // Blocul I: pentru variabila "x1" x1 = ls1; increment = fabs(ld1-ls1)/n; while(x1


68

printf("\n REZULTAT FINAL:\n\n"); fprintf(fp, "\n REZULTAT FINAL:\n\n"); printf(" x1Min=%lf\n x2Min=%lf\n\n fMin=%lf\n\n", x1Min, x2Min, fMin); fprintf(fp, " x1Min=%lf\n x2Min=%lf\n\n fMin=%lf\n\n", x1Min, x2Min, fMin); fclose(fp); }

Dei este o metod de optimizare de ordinul zero, pentru c nu folosete n calcule derivata de ordinul nti a funciei obiectiv, este aproape la fel de eficient ca i o metod de ordinul unu. Programul surs n limbaj C, dat mai sus, este construit pentru funcii dependente de dou variabile, dar trecerea la mai mult de dou variabile este foarte uoar. n plus, procedura de cutare efectiv este scris pe blocuri pentru a facilita implementarea mai multor variabile, n cazul n care este necesar. Totodat, trebuie avut n vedere definiia funciei obiectiv, precum i modul de apelare a acesteia. Pentru determinarea minimului funciei (1.1), intervalele de cutare au fost alese [-5,5] pentru variabila x1 i respectiv [-4,4] pentru variabila x2. Pe fiecare dintre cele dou axe s-au ales un numr de 1000 diviziuni, iar factorul de convergen a fost ales 1e-4. Dup rularea programului s-au obinut rezultatele din tabelul 3.4. Tabelul 3.4

iteraia x1Min x2Min fMin 1 4.960000 -2.368000 26.965334 2 2.320000 -1.728000 5.751951 3 1.000000 -1.312000 0.486993 4 0.360000 -1.056000 -0.751306 5 0.040000 -0.896000 -1.035328 6 -0.080000 -0.832000 -1.081494 7 -0.160000 -0.800000 -1.094400 8 -0.200000 -0.800000 -1.094400

Dup cum se poate observa din tabelul 3.4, dup 8 iteraii, s-au obinut pentru punctul de optim coordonatele:

;80.020.0min

X

iar funcia obiectiv are n acest punct valoarea: 0944.1min =f .


69

3.4 Metode de ordinul nti Din punct de vedere al numrului de iteraii necesare, sau al timpului maxim necesar atingerii soluiei cu o anumit precizie impus iniial, metodele de ordinul nti sunt mai eficiente dect cele de ordinul zero. Metodele de ordinul nti utilizeaz derivata de ordinul nti a funciei pentru determinarea vectorului direciei de minimizare. 3.4.1 Metoda pailor descendeni Una dintre cele mai cunoscute metode de optimizare, dar n acelai timp avnd i cele mai mici performane dintre metodele de ordinul nti, este metoda pailor descendeni, sau metoda gradientului. Totui, aceast metod este foarte important pentru c constituie punctul de plecare pentru alte metode de ordinul nti, mai sofisticate. Ca i n cadrul altor metode de optimizare se pleac dintr-un punct de start X0, ales la ntmplare n spaiul de cutare. Vectorul direcie de minimizare are aceeai

f2 0 2

2

0

2

7

6

6

6

5

5

5

4

4

4

4

33

3

33

3

2

2

2

2

2

2

1

1

1 1

1

1

0

0

0

00

0

1 1

X0

X1

X2

X*

F

F

FS

S

S

X1

X2

O

Fig.3.2 Interpretarea geometric a metodei pailor descendeni.


70

direcie ca i gradientul, dar sensul este contrar acestuia (vezi fig.3.2). Prin urmare, la iteraia q avem: ( );qq XFS = Acest vector este folosit n ecuaia (3.9), pentru procedura de cutare unidimensional. Metoda are o convergen lent, explicabil prin faptul c la determinarea direciei de minimizare, ncepnd de la iteraia a doua, nu se utilizeaz nici o informaie referitoare la iteraiile precedente. Dar, direcia pailor descendeni este folosit ca direcie iniial n majoritatea metodelor mai performante [Vanderplaats,1989]. n continuare este ilustrat programul surs n limbaj C, pentru metoda pailor descendeni. Ca de obicei, s-a folosit drept funcie test, funcia (1.1). Programul a fost construit pentru o funcie dependent de dou variabile, dar cu un efort minim din partea unui utilizator, se poate trece la funcii de n variabile.

// macrodefinitia functiei obiectiv #define f(u,v) 2*u+v+(u*u-v*v)+(u-v*v)*(u-v*v) #include #include void main( void ) { // punctul de start double x1zero, x2zero; // coordonatele punctului curent double x1unu, x2unu; // coordonatele vectorului directie de minimizare double s1, s2; double x1Min, x2Min, fMin; double epsilon = 1.0e-6; // increment pentru calculul derivatelor partiale double delta = 0.00001; double f_actual, f_precedent, f, deltaF = 1e3; double fa, fb, x1, x2, diff; double x1p, x2p; // factorul de salt initial double alfa0 = 0.001, alfa; // incrementul lui alfa double incAlfa = 0.001; FILE *fp; // deschidem fisierul "date_out.txt"

(3.10)


71

fp = fopen( "date_out.txt", "w" ); printf("\n\n METODA PASILOR DESCENDENTI\n"); fprintf(fp, "\n\n METODA PASILOR DESCENDENTI\n"); printf(" ========================================\n\n"); fprintf(fp, " ========================================\n\n"); // calculul functiei in punctul de start x1zero = 3.0; x2zero = 3.5; printf(" Coordonatele punctului de start:\n\n"); fprintf(fp, " Coordonatele punctului de start:\n\n"); printf(" x1zero=%lf x2zero=%lf\n\n", x1zero, x2zero); fprintf(fp, " x1zero=%lf x2zero=%lf\n\n", x1zero, x2zero); f_precedent = f(x1zero,x2zero); printf(" F(X0)=%lf\n\n", f_precedent); fprintf(fp, " F(X0)=%lf\n\n", f_precedent); fMin = f_precedent; // Calculul directiei de minimizare in X0 x1 = x1zero; x2 = x2zero; fa = f(x1,x2); x1 = x1 + delta; fb = f(x1,x2); diff = fa-fb; s1 = -diff/(-delta); x2 = x2 + delta; fb = f(x1,x2); diff = fa-fb; s2 = -diff/(-delta); // calculul lui X1 alfa = alfa0; x1unu = x1zero + alfa*s1; x2unu = x2zero + alfa*s2; f_actual = f(x1unu,x2unu); while(f_actual


72

fMin = f_actual; x1Min = x1unu; x2Min = x2unu; printf("s1=%lf s2=%lf x1Min=%lf x2Min=%lf F(X1)=%lf\n", s1, s2, x1Min, x2Min, f_actual); deltaF = fabs(f_precedent-f_actual); while(deltaF>=epsilon) { f = f_actual; alfa = alfa0; x1p = x1Min; x2p = x2Min; // Calculul directiei de minimizare in punctul actual x1 = x1p; x2 = x2p; fa = f(x1,x2); x1 = x1 + delta; fb = f(x1,x2); diff = fa-fb; s1 = -diff/(-delta); x2 = x2 + delta; fb = f(x1,x2); diff = fa-fb; s2 = -diff/(-delta); // calculul lui X1 alfa = alfa0; x1unu = x1p + alfa*s1; x2unu = x2p + alfa*s2; f_actual = f(x1unu,x2unu); while(f_actual


73

deltaF = fabs(f-f_actual); printf("s1=%lf s2=%lf x1Min=%lf x2Min=%lf fMin=%lf\n", s1, s2, x1Min, x2Min, f_actual); fprintf(fp, "s1=%lf s2=%lf x1Min=%lf x2Min=%lf fMin=%lf\n", s1, s2, x1Min, x2Min, f_actual); } printf("\n REZULTAT FINAL:\n\n"); fprintf(fp, "\n REZULTAT FINAL:\n\n"); printf(" x1Min=%lf\n x2Min=%lf\n\n fMin=%lf\n\n", x1Min, x2Min, fMin); fprintf(fp, " x1Min=%lf\n x2Min=%lf\n\n fMin=%lf\n\n", x1Min, x2Min, fMin); fclose(fp); }

La fel ca i-n cazul precedent, pentru determinarea minimului funciei (1.1), intervalele de cutare au fost alese [-5,5] pentru variabila x1 i respectiv [-4,4] pentru variabila x2, iar factorul de convergen a fost ales 1e-4. Dup rularea programului s-au obinut rezultatele din tabelul 3.5. Tabelul 3.5

Direcia de minimizare Minimul la iteraia i i s1 s2 x1 x2 Fmin

1 -8.114028 -6.229671 3.149386 1.798762 14.788042 2 -8.126472 -6.148920 3.141259 1.792613 14.734392 3 -8.138132 -6.070682 -0.260480 -0.744932 -1.088083 4 0.151748 0.091508 -0.260328 -0.744841 -1.088100 5 0.150868 0.089654 -0.260178 -0.744751 -1.088117 6 0.149998 0.087823 -0.260028 -0.744663 -1.088134 7 0.149136 0.086015 -0.259878 -0.744577 -1.088151 8 0.148283 0.084228 -0.259730 -0.744493 -1.088167 9 0.147439 0.082463 -0.259583 -0.744410 -1.088184

10 0.146604 0.080720 -0.259436 -0.744330 -1.088200 11 0.145777 0.078997 -0.259290 -0.744251 -1.088216 12 0.144959 0.077296 -0.259145 -0.744173 -1.088231 13 0.144149 0.075615 -0.259001 -0.744098 -1.088247 14 0.143348 0.073955 -0.258858 -0.744024 -1.088262 15 0.142554 0.072315 -0.258715 -0.743951 -1.088277 16 0.141769 0.070695 -0.258574 -0.743881 -1.088292 17 0.140991 0.069094 -0.258433 -0.743812 -1.088307 18 0.140222 0.067513 -0.258292 -0.743744 -1.088322


74

19 0.139460 0.065951 -0.258153 -0.743678 -1.088336 20 0.138706 0.064408 -0.258014 -0.743614 -1.088351 21 0.137959 0.062884 -0.257876 -0.743551 -1.088365 22 0.137221 0.061379 -0.257739 -0.743490 -1.088379 23 0.136489 0.059891 -0.257602 -0.743430 -1.088393 24 0.135765 0.058422 -0.257467 -0.743371 -1.088407 25 0.135048 0.056971 -0.257332 -0.743314 -1.088421 26 0.134339 0.055537 -0.257197 -0.743259 -1.088434 27 0.133636 0.054121 -0.257064 -0.743205 -1.088448 28 0.132941 0.052722 -0.256931 -0.743152 -1.088461 29 0.132252 0.051341 -0.256798 -0.743101 -1.088474 30 0.131571 0.049976 -0.256667 -0.743051 -1.088488 31 0.130896 0.048627 -0.256536 -0.743002 -1.088501 32 0.130228 0.047295 -0.256406 -0.742955 -1.088514 33 0.129566 0.045979 -0.256276 -0.742909 -1.088527 34 0.128263 0.043396 -0.256019 -0.742821 -1.088552 35 0.127621 0.042128 -0.255891 -0.742778 -1.088565 36 0.126985 0.040875 -0.240399 -0.737792 -1.089325 37 0.050250 -0.109504 -0.218943 -0.784550 -1.093692 38 0.106788 0.156332 -0.212001 -0.774389 -1.094302 39 0.047340 0.012721 -0.205705 -0.772697 -1.094420 40 0.016920 -0.047359 -0.199986 -0.788704 -1.094986 41 0.044032 0.059945 -0.197212 -0.784927 -1.095076 42 0.021050 0.004731 -0.194139 -0.784237 -1.095102 43 0.006589 -0.023636 -0.192261 -0.790973 -1.095210 44 0.020300 0.026034 -0.191002 -0.789359 -1.095227 45 0.010164 0.001818 -0.189356 -0.789064 -1.095234 46 0.002648 -0.012732 -0.188707 -0.792184 -1.095259 47 0.009918 0.012004 -0.188092 -0.791439 -1.095262 48 0.005101 0.000568 -0.187118 -0.791331 -1.095265 49 0.000860 -0.007436 -0.186945 -0.792825 -1.095271 50 0.004904 0.005441 -0.186641 -0.792448 -1.095272

Dup efectuarea a 50 iteraii, valoarea minim gsit este fmin = -1.095272, iar coordonatele acestui minim sunt [-0.186641,-0.792488]. 3.4.2 Metoda direciei conjugate Metoda direciei conjugate, sau de gradient conjugat, necesit doar o simpl modificare a algoritmului pailor descendeni, avnd n schimb, o cretere spectaculoas a ratei de convergen. Este la fel de simplu de programat ca i metoda pailor descendeni. Direcia conjugat se calculeaz n modul urmtor:


75

( ) ;1+= qqqq SXFS iar scalarul q este dat de relaia: ( ) ( ) ;/ 212 = qqq XFXF unde q este indice iterativ. n continuare este ilustrat programul scris n limbaj C, al metodei gradientului conjugat. Acesta este derivat din programul anterior, al metodei pailor descendeni, n care dup prima iteraie, s-a modificat calculul direciei de minimizare n conformitate cu relaiile (3.11) respectiv (3.12). Programul este construit tot pentru funcia obiectiv (1.1), dependent de dou variabile de decizie, dar cu aceleai simple modificri, poate fi adaptat la funcii dependente de n variabile.

// macrodefinitia functiei obiectiv #define f(u,v) 2*u+v+(u*u-v*v)+(u-v*v)*(u-v*v) #include #include void main( void ) { // punctul de start double x1zero, x2zero; // coordonatele punctului curent double x1unu, x2unu; // coordonatele vectorului directie de minimizare double s1, s2; double s1p, s2p; double beta; double x1Min, x2Min, fMin; double epsilon = 1.0e-4; // increment pentru calculul derivatelor partiale double delta = 0.00001; double f_actual, f_precedent, f, deltaF = 1e3; double fa, fb, x1, x2, diff; double x1p, x2p; // factorul de salt initial double alfa0 = 0.001, alfa; // incrementul lui alfa double incAlfa = 0.001; FILE *fp;

(3.11)

(3.12)


76

// deschidem fisierul "date_out.txt" fp = fopen( "date_out.txt", "w" ); printf("\n\n METODA DIRECTIEI CONJUGATE\n"); fprintf(fp, "\n\n METODA DIRECTIEI CONJUGATE\n"); printf(" ========================================\n\n"); fprintf(fp, " ========================================\n\n"); // calculul functiei in punctul de start x1zero = 3.0; x2zero = 3.5; printf(" Coordonatele punctului de start:\n\n"); fprintf(fp, " Coordonatele punctului de start:\n\n"); printf(" x1zero=%lf x2zero=%lf\n\n", x1zero, x2zero); fprintf(fp, " x1zero=%lf x2zero=%lf\n\n", x1zero, x2zero); f_precedent = f(x1zero,x2zero); printf(" F(X0)=%lf\n\n", f_precedent); fprintf(fp, " F(X0)=%lf\n\n", f_precedent); fMin = f_precedent; // Calculul directiei de minimizare in X0 x1 = x1zero; x2 = x2zero; fa = f(x1,x2); x1 = x1 + delta; fb = f(x1,x2); diff = fa-fb; s1 = -diff/(-delta); x2 = x2 + delta; fb = f(x1,x2); diff = fa-fb; s2 = -diff/(-delta); // calculul lui X1 alfa = alfa0; x1unu = x1zero + alfa*s1; x2unu = x2zero + alfa*s2; f_actual = f(x1unu,x2unu); while(f_actual


77

f_actual = f(x1unu,x2unu); } fMin = f_actual; x1Min = x1unu; x2Min = x2unu; printf("s1=%lf s2=%lf x1Min=%lf x2Min=%lf F(X1)=%lf\n", s1, s2, x1Min, x2Min, f_actual); deltaF = fabs(f_precedent-f_actual); // directia din punctul de start = directie precedenta s1p = s1; s2p = s2; while(deltaF>=epsilon) { f = f_actual; alfa = alfa0; x1p = x1Min; x2p = x2Min; // calculul lui beta beta = ((sqrt(s1*s1+s2*s2))*(sqrt(s1*s1+s2*s2)))/((sqrt(s1p*s1p+s2p*s2p))*(sqrt(s1p*s1p+s2p*s2p))); // Calculul directiei de minimizare in punctul actual x1 = x1p; x2 = x2p; fa = f(x1,x2); x1 = x1 + delta; fb = f(x1,x2); diff = fa-fb; s1 = -diff/(-delta); x2 = x2 + delta; fb = f(x1,x2); diff = fa-fb; s2 = -diff/(-delta); // calculul lui X1 alfa = alfa0; x1unu = x1p + alfa*(s1+beta*s1p); x2unu = x2p + alfa*(s2+beta*s2p); f_actual = f(x1unu,x2unu); while(f_actual


78

f_precedent = f_actual; alfa = alfa + incAlfa; x1unu = x1p + alfa*(s1+beta*s1p); x2unu = x2p + alfa*(s2+beta*s2p); f_actual = f(x1unu,x2unu); } fMin = f_actual; x1Min = x1unu; x2Min = x2unu; deltaF = fabs(f-f_actual); s1p = s1; s2p = s2; printf("s1=%lf s2=%lf x1Min=%lf x2Min=%lf fMin=%lf\n", s1, s2, x1Min, x2Min, f_actual); fprintf(fp, "s1=%lf s2=%lf x1Min=%lf x2Min=%lf fMin=%lf\n", s1, s2, x1Min, x2Min, f_actual); } printf("\n REZULTAT FINAL:\n\n"); fprintf(fp, "\n REZULTAT FINAL:\n\n"); printf(" x1Min=%lf\n x2Min=%lf\n\n fMin=%lf\n\n", x1Min, x2Min, fMin); fprintf(fp, " x1Min=%lf\n x2Min=%lf\n\n fMin=%lf\n\n", x1Min, x2Min, fMin); fclose(fp); }

La fel ca i-n cazul metodei pailor descendeni, pentru determinarea minimului funciei (1.1), intervalele de cutare au fost alese [-5,5] pentru variabila x1 i respectiv [-4,4] pentru variabila x2, iar factorul de convergen a fost ales 1e-4. Dup rularea programului s-au obinut rezultatele din tabelul 3.6. Tabelul 3.6

Direcia de minimizare Minimul la iteraia i i s1 s2 x1 x2 Fmin

1 -8.114028 -6.229671 3.159886 1.685762 15.249801 2 -8.955978 -4.439566 3.142816 1.675092 15.145568 3 -8.959414 -4.352029 3.124900 1.666301 15.025887 4 -8.946504 -4.292160 3.106994 1.657657 14.906236 5 -8.932347 -4.235545 3.089115 1.649129 14.786891 6 -8.917230 -4.181655 -0.623596 -0.101649 -0.568436 7 0.515031 -0.430522 -0.631999 -0.106261 -0.568305


79

8 0.550557 -0.388546 -0.147156 -0.478937 -0.839195 9 -0.952634 -2.189185 -0.195808 -0.790842 -1.095115

10 0.034073 0.050207 -0.196726 -0.792981 -1.095026 Dup efectuarea a doar 10 iteraii, valoarea minim gsit este fmin = -1.095026, iar coordonatele acestui minim sunt [-0.196726,-0.792981]. 3.4.3 Alte metode de ordinul nti n vederea accelerrii convergenei, o serie de metoode de optimizare de ordinul nti, au la baz ideea crerii unui ir care, n decursul procesului de optimizare, aproximeaz inversa matricei hessiene. Aceasta pentru c, fa de metoda direciei conjugate a lui Fletcher i Reeves, unde se ine cont de istoria iteraiilor doar prin intermediul parametrului scalar , aceste metode memoreaz informaiile referitoare la iteraiile precedente, stocndu-le ntr-un ir n-dimensional. Aceste metode poart numele de metode de metric variabil. Direcia de cutare este definit sub forma: ( );qq XFHS = dar, trebuie fcut observaia c H nu este matricea hessian (matricea derivatelor pariale de ordinul al doilea), a lui F(X), ci doar o aproximare a acesteia, n decursul optimizrii. Datorit faptului c aceste metode genereaz o aproximare a matricei H-1, ele prezint o convergen de o caracteristic similar metodelor de ordinul al doilea. De aceea, metodele de variabil metric poart numele de cvasi-metode Newton. Exist mai multe variante ale acestor metode, dintre care cele mai cunoscute sunt DFP (Davidon-Fletcher-Powell), respectiv BFGS (Broydon-Fletcher-Goldfarb-Shanno). 3.5 Metode de ordinul al doilea Metoda clasic de optimizare, de ordinul al doilea, este metoda lui Newton. Se pleac de la dezvoltarea funciei obiectiv, n serie Taylor, pn la termenii de ordinul al doilea inclusiv ( vezi relaia 3.3) i care pentru iteraia q, are expresia:

( ) ( ) ( ) ( ) ;21 XXHXXXFXFXF qTqq ++

n care: ;1 qq XXX = +

(3.13)

(3.14)

(3.15)


80

Dac difereniem ecuaia (3.14) n raport cu X i impunem condiia de anulare a gradientului funciei obiectiv, vom obine: ( )[ ] ( );1 qq XFXHX = innd cont de (3.15), relaia de recuren de ordinul doi, devine: ( )[ ] ( );11 qqqq XFXHXX = + Dac comparm ultimul termen din aceast ecuaie, cu vectorul direcie S, din ecuaia (3.7), inversa matricei hessiene ( )[ ] 1qXH joac rolul factorului de salt optim, notat *. Principala dificultate a metodei lui Newton, const n calculul matricei hessiene, calcul ce poate fi laborios, situaie n care viteza de calcul va avea de suferit. De asemenea, exist posibilitatea ca matricea hessian s fie singular, sau s nu fie pozitiv definit, condiii ce garanteaz soluia de minim pentru funcia obiectiv. Dac funcia obiectiv este liniar n funcie de una sau mai multe variabile de decizie, ntotdeauna matricea hessian va fi singular, iar direcia S poate s fie nedeterminat, ceea ce va conduce la un rezultat eronat. Uneori, datorit pailor mari dai de metoda lui Newton, putem ajunge la situaia unor oscilaii n jurul soluiei, aa numitul efect Maratos [Marusciac,1973]. De aceea este necesar la fiecare iteraie determinarea unor valori rezonabile pentru mrimea acestor pai. Dar dac matricea hessian se poate calcula uor, atunci ntotdeauna, metoda lui Newton este cea mai recomandat. 3.6 ALGORITMUL DE CONVERGEN

Un aspect important din cadrul oricrui algoritm de optimizare l reprezint detectarea momentului n care ne-am apropiat de soluia optim cu precizia , cu alte cuvinte determinarea convergenei. De regul, se folosesc mai multe criterii de verificare a convergenei. Un prim criteriu este reprezentat de diferena absolut a valorilor funciei obiectiv. Pe msur ce procesul de cutare progreseaz n direcia punctului de optim, diferena absolut dintre valorile funciei obiectiv ntre dou iteraii succesive devine tot mai mic (vezi fig.3.3). Dac notm diferena absolut a valorilor funciei obiectiv cu dF1, atunci la iteraia q: ( ) ( ) ;11 = qq XFXFdF Vom considera ndeplinit condiia de convergen atunci cnd dF1 a, unde a reprezint o cantitate foarte mic (echivalent cu precizia de calcul impus iniial).

(3.16)

(3.17)

(3.18)


81

Deoarece, se poate ntmpla ca pe anumite poriuni, funcia obiectiv s prezinte o pant lent, iar apoi panta s creasc din nou (vezi fig.3.4), exist riscul ca odat ndeplinit condiia (3.18), s deducem c punctul cutat este xi+1i nu x*. De aceea condiia de convergen trebuie ndeplinit de mai multe ori, consecutiv, pentru a ne asigura c am determinat adevratul punct de optim. De asemenea, se poate ntmpla ca n zona punctului de optim panta funciei obiectiv s fie foarte mic pe poriuni mari, astfel nct procesul de cutare s se opreasc prea devreme. Pentru a evita acest lucru, putem s atribuim factorului de convergen a o valoare foarte mic. Totui nu se recomand aceasta, deoarece o valoare foarte mic a factorului de convergen chiar dac conduce la creterea preciziei de calcul, duce la creterea numrului de iteraii i deci implicit a numrului calculelor, respectiv a timpului de lucru pe calculator.

Pentru a nu lucra cu valori foarte mici ale lui a, putem folosi drept factor de convergen diferena relativ a valorilor funciei obiectiv ntre dou iteraii succesive. Astfel:

;}10),(max{ 10

12 = qXF

dFdF

n relaia (3.19), drept numitor folosim valoarea maxim dintre valoarea funciei

Fi

Fj

F

xi xi+1 xj xj+1 x

Fig.3.3 Variaia diferenei absolute a valorilor funciei obiectiv, ntre dou iteraii succesive.

(3.19)


82

obiectiv la iteraia cu numrul q i respectiv o cantitate foarte mic (de ex. 10-10), pentru a evita mprirea cu zero pe calculator. La fel ca i n cazul precedent, vom considera ndeplinit condiia de convergen dac dF2 r, unde r este o cantitate mica, impus iniial. De asemenea, aceast condiie trebuie ndeplinit de mai multe ori consecutiv, din aceleai motive, precizate mai sus.

n cazul n care putem determina derivata de ordinul nti a funciei obiectiv (sau gradientul funciei obiectiv, dac funcia depinde de mai multe variabile de decizie), pe cale numeric (cu diferene finite), atunci valoarea derivatei (sau gradientului) devine la rndul ei un criteriu de convergen. Din punct de vedere geometric, derivata de ordinul nti a funciei obiectiv ntr-un anumit punct, este reprezentat de tangenta la graficul funciei n acel punct. Astfel, n punctul de optim (minim sau maxim), tangenta la graficul funciei obiectiv este n poziie orizontal. Prin urmare, este suficient s verificm momentul n care valoarea derivatei de ordinul nti (sau valoarea gradientului funciei obiectiv) se apropie de zero cu o anumit precizie, adic:

( ) ;gradqXF

F

xi xi+1 x* x

F(xi) F(xi+1)

O

(3.20)

Fig.3.4 Variaii posibile ale pantei funciei obiectiv, traversate n decursul procesului de cutare a punctului de optim.


83

Un criteriu des folosit n cadrul algoritmilor de convergen este numrul maxim de iteraii, notat qmax. Dei foarte important, acest criteriu nu ne furnizeaz informaii legate de poziia curent din procesul de cutare n raport cu poziia punctului de optim. El este util atunci cnd algoritmul ales pentru determinarea soluiei optime are o convergen mult prea lent, sau programul de optimizare conine o greeal de programare, ceea ce conduce la repetarea unor calcule la nesfrit. Dm n cele ce urmeaz, n pseudocod, algoritmul de verificare a convergenei.

// intrare din programul de optimizare q = 0; n1=0, n2=0, n3=0; N=5; while( q


84

n3 = 0; } // revenire n programul de optimizare

Optimizare fara restrictii

Documents

Transcript of Optimizare fara restrictii