LabSem

Metode de Optimizare NumericaCulegere de probleme

Ion Necoara Dragos Clipici Andrei Patrascu

Departamentul de Automatica si Ingineria SistemelorUniversitatea Politehnica din Bucuresti

Email: {ion.necoara,dragos.clipici,andrei.patrascu}@acse.pub.ro

2013

Prefata

Lucrarea de fata este construita dupa structura cursului de Tehnicide optimizare, predat de primul autor la Facultatea de Automatica siCalculatoare a Universitatii Politehnica din Bucuresti si se adreseazastudentilor de la facultatile cu profil tehnic sau cu tematici adiacenteingineriei sistemelor si calculatoarelor. Aceasta contine o serie deprobleme rezolvate de laborator si seminar, avand rolul de suportpentru aplicatiile teoretice/practice abordate la diferite materii ce includnotiuni din teoria optimizarii. La nceputul fiecarui capitol este introdusun succint breviar teoretic n care se definesc instrumentele necesarentregului capitol, iar la final este expusa o sectiune cu problemepropuse ce au ca scop posibilele extensii ale aplicatiilor rezolvate.Rezolvarea problemelor de laborator a fost implementata n mediul deprogramare Matlab, iar n corpul rezolvarii au fost incluse secventade cod corespunzatoare si figuri ce reflecta performantele practice alealgoritmilor folositi. Pe tot parcursul lucrarii prezentam multe exemplede aplicatii si exercitii rezolvate pentru a face mai accesibila ntelegereateoriei si conceptelor legate de optimizarea numerica. Precizam ca oparte din algoritmii prezentati nca reprezinta o tema de cercetare deactualitate n domeniul optimizarii n care autorii activeaza. Capitoleleau fost structurate gradat, dupa dificultatea si importanta notiunilorteoretice necesare n cadrul aplicatiilor, n numar de sapte:

Cap 1: Functii ale pachetului de optimizare Matlab;Cap 2: Optimizare convexa;Cap 3: Algoritmi de ordinul I pentru probleme neconstranse;Cap 4: Algoritmi de ordinul II pentru probleme neconstranse;Cap 5: Optimizare constransa;Cap 6: Algoritmi pentru probleme constranse;Cap 7: Aplicatii din inginerie.

Desi fiecare dintre capitole dezvolta cazuri concrete de aplicatii alealgoritmilor n cauza, ultimul este alocat unei serii de aplicatii reale, deactualitate, din diverse domenii ale ingineriei, precum reglarea (controlul)optimala, prelucrarea semnalelor sau nvatarea automata. Fiecareaplicatie ilustreaza tehnici specifice de modelare matematica a uneiprobleme si strategii de formulare a modelului matematic de optimizare

6pentru aplicatia respectiva. Acest capitol final demonstreaza astfelaplicabilitatea metodelor numerice de optimizare prezentate n capitoleleanterioare n probleme extrem de studiate din inginerie.1In final, precizam ca am ncercat sa reducem din caracterul tehnic allimbajului pentru o expunere cat mai clara si n acelasi timp, pentruca cititorul sa si poata fixa cat mai usor notiunile fundamentale deoptimizare.

Autorii, mai 2013

Cuprins

1 Functii ale pachetului de optimizare MATLAB 11

1.1 Preliminarii . . . . . . . . . . . . . . . . . . . . . . . . . 111.1.1 Mediul Matlab . . . . . . . . . . . . . . . . . . . 12

1.2 Probleme rezolvate de laborator . . . . . . . . . . . . . . 131.2.1 Minimizare neconstransa . . . . . . . . . . . . . . 131.2.2 Minimizare constransa . . . . . . . . . . . . . . . 171.2.3 Programare liniara . . . . . . . . . . . . . . . . . 241.2.4 Programare patratica . . . . . . . . . . . . . . . . 27

1.3 Probleme rezolvate de seminar . . . . . . . . . . . . . . . 291.3.1 Diferentiere multivariabila . . . . . . . . . . . . . 291.3.2 Probleme de optimizare generale . . . . . . . . . 30

1.4 Probleme propuse . . . . . . . . . . . . . . . . . . . . . . 35

2 Probleme de optimizare convexa 38

2.1 Preliminarii . . . . . . . . . . . . . . . . . . . . . . . . . 382.1.1 Multimi convexe . . . . . . . . . . . . . . . . . . 382.1.2 Functii convexe . . . . . . . . . . . . . . . . . . . 40

2.2 Pachet CVX . . . . . . . . . . . . . . . . . . . . . . . . . 412.2.1 Instalare . . . . . . . . . . . . . . . . . . . . . . . 422.2.2 Elemente de baza . . . . . . . . . . . . . . . . . . 43

2.3 Probleme rezolvate de laborator . . . . . . . . . . . . . . 442.4 Probleme rezolvate de seminar . . . . . . . . . . . . . . . 48

2.4.1 Multimi si functii convexe . . . . . . . . . . . . . 482.5 Probleme propuse . . . . . . . . . . . . . . . . . . . . . . 57

3 Metode de ordinul I 63

3.1 Preliminarii . . . . . . . . . . . . . . . . . . . . . . . . . 633.2 Probleme rezolvate de laborator . . . . . . . . . . . . . . 64

3.2.1 Metoda Gradient . . . . . . . . . . . . . . . . . . 643.2.2 Metoda gradientilor conjugati . . . . . . . . . . . 70

8 Cuprins

3.3 Probleme rezolvate de seminar . . . . . . . . . . . . . . . 723.3.1 Aproximari patratice . . . . . . . . . . . . . . . . 723.3.2 Conditii de optimalitate de ordin I . . . . . . . . 733.3.3 Metoda gradient . . . . . . . . . . . . . . . . . . 75


4 Metode de ordinul II 85


4.2.1 Metoda Newton . . . . . . . . . . . . . . . . . . . 874.2.2 Metode cvasi-Newton . . . . . . . . . . . . . . . . 90

4.3 Probleme rezolvate de seminar . . . . . . . . . . . . . . . 924.3.1 Metoda Newton si metoda BFGS . . . . . . . . . 92


5 Probleme de optimizare constransa 103


5.2.1 Formularea unei probleme de optimizare n formastandard . . . . . . . . . . . . . . . . . . . . . . . 104

5.2.2 Calcularea proiectiei ortogonale a unui punct pe omultime convexa . . . . . . . . . . . . . . . . . . 106

5.2.3 Metoda Gauss-Newton . . . . . . . . . . . . . . . 1115.2.4 Metoda gradientului proiectat . . . . . . . . . . . 113

5.3 Probleme rezolvate de seminar . . . . . . . . . . . . . . . 1155.3.1 Problema duala, conditii si puncte Karush-Kuhn-

Tucker . . . . . . . . . . . . . . . . . . . . . . . . 1155.4 Probleme propuse . . . . . . . . . . . . . . . . . . . . . . 129

6 Metode pentru probleme de optimizare constransa 132

6.1 Preliminarii . . . . . . . . . . . . . . . . . . . . . . . . . 1326.1.1 Probleme neconvexe cu constrangeri de egalitate . 1326.1.2 Probleme convexe cu constrangeri de egalitate . . 1336.1.3 Probleme convexe generale . . . . . . . . . . . . . 134

6.2 Probleme rezolvate de laborator . . . . . . . . . . . . . . 1376.2.1 Metoda Lagrange-Newton . . . . . . . . . . . . . 1376.2.2 Metoda Newton extinsa . . . . . . . . . . . . . . 1386.2.3 Metoda de punct interior . . . . . . . . . . . . . . 140

6.3 Probleme rezolvate de seminar . . . . . . . . . . . . . . . 145

Cuprins 9

6.3.1 Metoda Newton pentru probleme constranse . . . 1456.3.2 Probleme de control optimal . . . . . . . . . . . . 147


7 Aplicatii din inginerie 157

7.1 Control optimal . . . . . . . . . . . . . . . . . . . . . . . 1577.1.1 Control optimal aplicat unui robot E-Puck . . . . 1587.1.2 Control optimal aplicat unei instalatii cu patru

rezervoare . . . . . . . . . . . . . . . . . . . . . . 1647.2 Problema Google . . . . . . . . . . . . . . . . . . . . . . 1677.3 Clasificarea de imagini . . . . . . . . . . . . . . . . . . . 173

Bibliografie 179

10 Cuprins

Capitolul 1

Functii ale pachetului de

optimizare MATLAB

1.1 Preliminarii

Explozia informationala din ultimele patru decenii a condus la algoritminumerici de optimizare indispensabili pentru abordarea diferitelorprobleme din domeniile ingineriei, fizicii, analizei numerice, teorieisistemelor si controlului, prelucrarii semnalelor, nvatarii automateetc. Precizand ca aceasta scurta enumerare reprezinta doar cele maievidente domenii de aplicatie ale algoritmilor de optimizare, constatam cateoria optimizarii este considerata un instrument valoros n majoritateadomeniilor stiintifice.In prezent exista algoritmi eficienti (rapizi) destinati rezolvariiproblemelor de optimizare convexe/neconvexe cu dimensiuni mici si mediintr-un interval de timp relativ scurt. Insa multe dificultati si ntrebarifara raspuns apar n cazul problemelor cu dimensiuni mari (e.g. 107109variabile). Dimensiunile problemelor ntalnite n domeniile enumerateanterior se afla n crestere continua, astfel ncat se formuleaza relativ desprobleme cu dimensiuni de zeci de milioane de variabile. Ca urmare,algoritmii de optimizare raman un subiect de actualitate n cercetareamoderna din matematica aplicata.In continuare vom prezenta si formula principalele clase de probleme deoptimizare si n plus, functiile Matlab uzuale pentru rezolvarea numericaa acestor tipuri de probleme.

12 Capitolul 1. Functii ale pachetului de optimizare MATLAB

1.1.1 Mediul Matlab

Matlab (Matrix laboratory) este un mediu de calcul numeric si unlimbaj de programare de generatia a patra, dezvoltat de companiaMathworks. Mediul Matlab ofera o interfata usor de utilizat pentruoperatii cu matrice, implementare de algoritmi, grafice de functii etc.

Figura 1.1: Matlab.

Inca de la sfarsitul anilor 1970, Matlab a castigat o audienta larga ncadrul academic si n comunitatea matematicii aplicate. Incepand cudezvoltarea de pachete de programe destinate problemelor numerice dealgebra liniara, au fost incluse si dezvoltate ulterior functii (toolbox-uri)specializate n rezolvarea de probleme de optimizare, rezolvarea desisteme de ecuatii, modelare si identificare de sisteme, control etc.Pachetul Optimization Toolbox, din cadrul mediului Matlab, furnizeazaalgoritmi utilizati la scara larga, destinati problemelor de optimizarestandard si de mari dimensiuni. De asemenea, pachetul include si functiice rezolva probleme de optimizare structurate: programare liniara,programare patratica, programare binara, probleme CMMP neliniare,rezolvare de sisteme de ecuatii neliniare si optimizare multicriteriala.Dificultatea problemelor de optimizare ce apar de cele mai multe orin practica limiteaza capacitatea algoritmilor existenti la gasirea uneisolutii aproximative corespunzatoare problemei date. De exemplu,

1.2. Probleme rezolvate de laborator 13

problema de optimizare ce presupune gasirea punctului global de maximcorespunzator functiei sinc(x) = sinx

xeste o problema dificil de rezolvat.

Pentru motivarea acestei afirmatii se poate observa usor n Fig. 1.2 cafunctia sinc(x) prezinta o infinitate de maxime locale (denivelarile dela baza dealului), despre care nu cunoastem a priori daca sunt sau nuvalori maxime globale. Pentru a ilustra simplitatea cu care putem figurafunctiile n mediul Matlab, prezentam n continuare secventa de cod caregenereaza graficul din Fig. 1.2.

[X,Y] = meshgrid(-10:0.25:10,-10:0.25:10);

f = sinc(sqrt((X/pi).^2+(Y/pi).^2));

surf(X,Y,f);

axis([-10 10 -10 10 -0.3 1])

xlabel({\bfx})

ylabel({\bfy})

zlabel({\bfsinc} ({\bfR}))

Figura 1.2: Graficul functiei sinc(x) .

1.2 Probleme rezolvate de laborator

1.2.1 Minimizare neconstransa

Problemele de optimizare neconstransa presupun minimizarea (saumaximizarea) unei anumite functii obiectiv (criteriu) prin intermediulunui set de variabile de decizie, fara ca acestea sa respecte anumiterestrictii (constrangeri). Se considera functia (n general, neliniara)


f : Rn R de doua ori diferentiabila, careia i se asociaza problemade optimizare:

minxRn

f(x),

unde x are rolul vectorului variabilelor de decizie. Problema deoptimizare se considera rezolvata daca s-a obtinut un vector x pentrucare valoarea functiei f() n x este minima, adica f(x) f(x)pentru orice x Rn. Multimea n care se efectueaza cautarea senumeste multime fezabila, care n cazul de fata este data de ntreg spatiulEuclidian Rn.

Functii MATLAB

Pentru cazul n care functia obiectiv nu are o structura prestabilitasau o forma particulara, sunt disponibile urmatoarele functii MATLABpentru aproximarea unui punct de minim local n cazul unei problemeneconstranse:

1. Functia [...]=fminunc(...) returneaza un punct de minim localal functiei introduse ca argument, cu ajutorul informatiei de ordinul I(gradient) si a celei de ordin II (Hessiana). Sintaxa functiei are expresia:

x=fminunc(fun,x0)

x=fminunc(fun,x0,optiuni)

[x,val]=fminunc(...)

[x,fval,exitflag,output,grad,hessian]=fminunc(...)

Variabila fun reprezinta functia obiectiv diferentiabila ce trebuiefurnizata ca o variabila de tip function handle. O variabila de tipfunction handle poate fi, e.g. un fisier objfun.m care primestela intrare variabila de decizie x si returneaza valoarea functiei nx. Daca metoda de rezolvare a problemei necesita si gradientulfunctiei, iar optiunea GradObj este setata on, atunci argumentulfun furnizeaza, n plus, si valoarea gradientului n punctul x. Dacametoda de rezolvare a problemei necesita Hessiana, iar optiuneaHessian este setata on, atunci argumentul fun furnizeaza, nplus, si valoarea hessianei n punctul x.

Vectorul x0 reprezinta punctul initial de unde porneste procesulde cautare al minimului si trebuie furnizat ca un vector din Rn.


Variabila optiuni reprezinta setul de optiuni specific fiecarei rutineMatlab. Comanda optiuni=optimset(fminunc) afiseaza setulde optiuni implicit. Fiecare dintre acestea se poate modifica nfunctie de problema de minimizare. De exemplu, pentru setareaoptiunii GradObj pe on, comanda este optiuni.GradObj=on.

Pentru o documentare amanuntita introduceti urmatoarea comanda:help fminunc

2. Functia [...]=fminsearch(...) returneaza un punct de minim localal functiei introduse ca argument, utilizand nsa o metoda diferita de ceafolosita de fminunc, anume metoda Nedler-Mead, ce utilizeaza informatiede ordin zero (evaluarea functiei). Sintaxa este similara celei prezentatepentru functia fminunc.

Exemplul 1. Consideram problema de optimizare

minxR2

f(x)(= ex1(4x21 + 2x

22 + 4x1x2 + 2x2 + 1)

).

Sa se rezolve problema (sa se gaseasca un punct de minim local) utilizandrutina fminunc, cu punctul de initializare x0 = [1;1].Rezolvare. Se creeaza fisierul objfun.m ce va contine urmatoareasecventa de cod:

function f=objfun(x)

f = exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);

Se apeleaza rutina fminunc n urmatorul mod:

x0=[-1,1]; optiuni=optimset(LargeScale,off);

[x,fval,exitflag,output]=fminunc(@objfun,x0,optiuni);

Exemplul 2. Consideram problema de optimizare:

minxR

f(x) (= cos 3x+ sin 2x) .

(i) Sa se rezolve problema (sa se gaseasca un punct de minim local)utilizand rutina fminunc, cu punctul de initializare x0 = 0.

(ii) Sa se rezolve problema cu ajutorul metodei sectiunii de aur.


Rezolvare. (i) Se creeaza fisierul objfun.m cu urmatoarea secventa decod:


f = cos(3*pi*x) + sin(2*pi*x);

Se apeleaza rutina fminunc:

x0=0; optiuni=optimset(LargeScale,off);

[x,fval,exitflag,output]=fminunc(@objfun,x0,optiuni);

(ii) Deoarece functia este periodica, de perioada 2, putem aplica metodasectiunii de aur, n care se alege lungimea intervalului initial egala cu 2,a carei implementare este expusa n cele ce urmeaza. Se creeaza fisierulobjfun.m si goldsection.m:


f = cos(3*pi*x) + sin(2*pi*x);

end

function [xst,fst]= goldsection(x,eps)

a=x-pi; b=x+pi; lambda=a+0.382*(b-a);

miu=a+0.618*(b-a);

while ((b-a)>=eps)

if (objfun(lambda)


neliniara f : R R cu structura necunoscuta, dorim aproximareaacesteia cu un polinom g : R R,

g(z) = anzn + + a1z + a0.

Datele cunoscute n legatura cu functia neliniara f reprezinta setul deperechi: (y1, f(y1)), . . . , (ym, f(ym)), unde scalarii yi R sunt dati.Notand vectorul coeficientilor polinomului cu x = [a0 . . . an]

T si b =[f(y1) . . . f(ym)]

T , problema se reduce la gasirea vectorului optim decoeficienti ce satisface relatia:

Ax =

1 y1 . . . y

n1

1 y2 . . . yn2

......

......

1 ym . . . ynm

a0a1...an

f(y1)f(y2)...

f(ym)

= b.

Echivalent, obtinem urmatoarea problema de optimizare faraconstrangeri:

minxRn+1

f(x)(= Ax b2) . (1.1)

Sa se rezolve problema (1.1), pe cazuri particulare pentru matricea A sivectorul b de diferite dimensiuni, utilizand rutinele fminunc si lsqlin.Sa se compare rezultatele obtinute. Ce se observa?

1.2.2 Minimizare constransa

Fie functia (n general neliniara) f : Rn R. Forma generala a uneiprobleme de optimizare cu constrangeri (numite si restrictii), asociatefunctiei f , se formuleaza n urmatorul mod:

minxRn

f(x)

s.l.: Cx d, g(x) 0, l x u, (1.2)Ax = b, h(x) = 0,

unde l, u Rn, d Rm1 , b Rp1, C Rm1n, A Rp1n, iar g(x) :Rn Rm2 , h(x) : Rn Rp2 sunt functii multidimensionale (vectori de

functii), reprezentand constrangerile neliniare de inegalitate si respectiv,de egalitate.


Exemplul 4. Consideram urmatoarea problema de optimizare:

minxR2

(x1 3)2 + (x2 2)2

s.l.: x21 x2 3 0, x2 1 0, x1 0.

Functia obiectiv si cele trei constrangeri de inegalitate sunt definite deurmatoarele expresii:

f(x1, x2) = (x1 3)2 + (x2 2)2,

g1(x1, x2) = x21 x2 3, g2(x1, x2) = x2 1, g3(x1, x2) = x1.

Fig 1.3 ilustreaza multimea fezabila. Problema se reduce la a gasi unpunct n multimea fezabila n care (x1 3)2 + (x2 2)2 ia cea mai micavaloare. Observam ca punctele [x1 x2]

T cu (x1 3)2+ (x2 2)2 = c suntcercuri de raza c cu centrul n [3 2]T . Aceste cercuri se numesc multimilenivel sau contururile functiei obiectiv avand valoarea c. Pentru aminimiza c trebuie sa gasim cercul cu cea mai mica raza care intersecteazamultimea fezabila. Dupa cum se observa din Fig. 1.3, cel mai mic cerccorespunde lui c = 2 si intersecteaza multimea fezabila n punctul deoptim x = [2 1]T .

5 4 3 2 1 0 1 2 3 4 55

4

3

2

1

0

1

2

3

4

5

g1

g3 contururi functie obiectiv

punct optim (2,1)zona fezabila

g2(3,2)

Figura 1.3: Solutia grafica a problemei de optimizare.


Functii MATLAB

Pentru cazul n care functia de minimizat (maximizat) nu are o structuraprestabilita sau o forma particulara, pentru rezolvarea problemei deoptimizare (1.2) este disponibila functia MATLAB fmincon ce prezintaurmatoarea sintaxa:

x = fmincon(fun,x0,C,d,A,b,l,u,nonlcon,optiuni)

[x,fval,exitflag,output,lambda,grad,hessian] = fmincon(...).

variabila fun reprezinta functia obiectiv ce trebuie furnizata cao variabila de tip function handle. Daca metoda de rezolvare aproblemei cere si gradientul functiei, iar optiunea GradObj estesetata on, atunci argumentul fun furnizeaza, n plus, si valoareagradientului n punctul x. Daca metoda de rezolvare a problemeicere Hessiana, iar optiunea Hessian este setata on, atunciargumentul fun furnizeaza, n plus, si valoarea Hessianei n punctulx;

matricele si vectorii C, d, A, b, l, u definesc constrangerile liniare custructura prezentata n modelul (1.2);

variabila de tip handle nonlcon returneaza constrangerile neliniarede egalitate si inegalitate reprezentate n (1.2);

x0 reprezinta punctul initial de unde porneste procesul de cautareal minimului si trebuie furnizat ca un vector din Rn;

variabila optiuni reprezinta setul de optiuni specific fiecarei rutineMATLAB. Comanda optiuni=optimset(fmincon) afiseaza setulde optiuni implicit. Fiecare dintre acestea se poate modifica nfunctie de problema de minimizare. De exemplu, pentru setareaoptiunii GradObj pe on, comanda este optiuni.GradObj=on.

Pentru o documentare amanuntita introduceti urmatoarea comanda:help fmincon

Exemplul 5. Consideram problema de optimizare:

minxR2

f(x)(= ex1(4x21 + 2x

22 + 4x1x2 + 2x2 + 1)

).

s.l.: g(x) 0({

x1x2 x1 x2 + 1.5 0x1x2 10 0

).


Sa se rezolve local problema utilizand rutina fmincon, cu punctul deinitializare x0 = [1 1]

T .

Rezolvare. Se creeaza fisierul objfun.m ce defineste functia obiectiv siconfun.m n care se definesc constrangerile:


f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);

function[g,h]=confun(x)

g=[x(1)*x(2)-x(1)-x(2)+1.5 -x(1)*x(2)-10];

h=[];

Apoi se apeleaza rutina pentru minimizare cu restrictii:

x0=[-1,1];

optiuni=optimset(fmincon);

optiuni.LargeScale=off;

[x,fval]=fmincon(@objfun,x0,[],[],[],[],...

[],[],@confun,optiuni)

Exemplul 6 (Analiza statistica). Analiza datelor si interpretareaacestora ntr-un sens cat mai corect este preocuparea principala dindomeniul statisticii. Problema se formuleaza n urmatorul mod: pe bazaunei colectii de date cunoscute (reprezentate n Fig. 1.4 prin puncte),sa se realizeze predictia cu o eroare cat mai mica a unui alt set de datepartial cunoscut. In termeni matematici, aceasta problema presupunedeterminarea unei directii de-a lungul careia elementele date (punctele)tind sa se alinieze, astfel ncat sa se poata prezice zona de aparitie apunctelor viitoare. S-a constatat ca directia de cautare este data devectorul singular al matricei formate din colectia de puncte date, ce poatefi gasit prin intermediul urmatoarei probleme de optimizare:

maxxRn

f(x)

(=

1

2xTATAx

)(1.3)

s.l.: g(x) 0 (x 1) ,unde A Rmn reprezinta matricea ale carei coloane [a1 . . . an] suntpunctele cunoscute initial. Considerand cazul particular n care se da

matricea A =

[1 0.5 00 2 1

], sa se rezolve problema precedenta cu ajutorul

functiei fmincon, alegand punctul initial x0 = [1 2 0.5]T .


2.5 2 1.5 1 0.5 0 0.5 1 1.5 2 2.510

5

0

5

10

15

Figura 1.4: Un exemplu de dispersie a datelor n spatiul R2. Observam capunctele se aliniaza de-a lungul dreptei ce iese n evidenta.

Rezolvare. Se creeaza fisierele objfun.m si confun.m pentru functiaobiectiv si respectiv, pentru constrangeri :

function [f]=objfun(x)

f=x(1)^2+4.25*x(2)^2+x(3)^2+0.5*x(1)*x(2)+2*x(2)*x(3);

function [g,h]=confun(x)

g=x(1)^2+x(2)^2 + x(3)^2 - 1;

h=[];

Se apeleaza rutina pentru minimizare cu restrictii:

x0=[-1,-2,0.5];



[x,fval]=fmincon(@objfun,x0,[],[],[],...

[],[],[],@confun,optiuni)

Exemplul 7 (Analiza stabilitatii robuste). Stabilitatea robusta asistemelor liniare afectate de perturbatii reprezinta o tema de cercetarede actualitate n teoria controlului. Numeroase aplicatii practice necesitaanaliza stabilitatii robuste (Fig. 1.5) pentru asigurarea functionariicorespunzatoare a diferitelor sisteme afectate de zgomot si perturbatii(e.g. sistemul de viraj al autovehiculelor, sistemul de ghidare al uneimacarale).


Figura 1.5: Exemplu de reglare robusta n cazul unui sistem cu parametriafectati de incertitudini.

In general, problema se reduce la studiul stabilitatii unui polinom cucoeficienti afectati de incertitudini. Consideram ca exemplu urmatorulpolinom:

s3 + (1 + a1 + a2)s2 + (a1 + a2 + 3)s+ (1 + 6a1 + 6a2 + 2a1a2), (1.4)

unde parametrii a1 si a2 satisfac conditiile:

a1 [1.4 1.1, 1.4 + 1.1] si a2 [0.85 0.85, 0.85 + 0.85],

unde > 0. Parametrii a1 si a2 pot reprezenta cantitati fizice (e.g.lungimi, mase) ce sunt afectate de incertitudini n mod obisnuit, sauparametrii regulatorului ce pot varia datorita erorilor de implementare.Suntem interesati n determinarea celei mai mari valori a variabilei pentru care polinomul ramane stabil. Aceasta problema se poatereformula ntr-una de optimizare folosind criteriul Hurwitz : un polinomde gradul 3 de forma

p(s) = s3 + q2s2 + q1s+ q0

este stabil daca si numai daca q2q1 q0 > 0 si q0 > 0. De aceea,consideram marginea maxima de stabilitate data de

= min{1, 2},


unde 1 si 2 sunt solutiile urmatoarelor probleme de optimizare:

1 = mina,

s.l.: 2 2a1 2a2 + a21 + a22 0,1.4 1.1 a1 1.4 + 1.10.85 0.85 a2 0.85 + 0.85,

2 = mina,

s.l.: 1 + 6a1 + 6a2 + 2a1a2 0,1.4 1.1 a1 1.4 + 1.10.85 0.85 a2 0.85 + 0.85.

Sa se rezolve problemele de optimizare neliniara precedente utilizandfunctia fmincon.

Rezolvare. Pentru a simplifica rezolvarea furnizam secventa de codMatlab ce rezolva prima problema de optimizare si determina 1,cea de-a doua rezolvandu-se n mod analog. Consideram urmatoareasecventa de cod aferenta functiei obiectiv si constrangerilor problemeide optimizare ce determina 1:

function [f]=objfun(delta)

f=delta;

end

function [g]=confun(x)

g=[1+6*x(1)+6*x(2)+2*x(1)*x(2);x(1)-1.4-1.1*x(3);...

-x(1)+1.4-1.1*x(3); x(2)-0.85-0.85*x(3);...

-x(2)+0.85-0.85*x(3)];

end

Se apeleaza apoi rutina pentru minimizare cu restrictii:

x0=[1.4 0.85 1000];



[x,fval]=fmincon(@objfun,x0,[],[],[],...

[],[],[],@confun,optiuni)


1.2.3 Programare liniara

Pentru cazurile particulare de probleme de optimizare (e.g. functieobiectiv liniara sau patratica), exista functii Matlab specifice ce rezolvaproblemele de optimizare aferente prin intermediul unor algoritmidedicati respectivelor clase de probleme.Functia linprog rezolva probleme de programare liniara, ce vizeazaminimizarea unei functii obiectiv liniare, n prezenta unui set deconstrangeri liniare. O problema de programare liniara se formuleazan general dupa cum urmeaza:

minxRn

cTx

s.l.: Cx d, l x u,Ax = b,

unde c, l, u Rn, d Rm, b Rp, C Rmn, A Rpn. Sintaxa functieilinprog cuprinde:

x = linprog(c,C,d,A,b,l,u,x0,optiuni)

[x,fval,exitflag,output,lambda] = linprog(. . . ).

Pentru mai multe detalii, se poate apela comanda: help linprog

Exemplul 8 (Repartizarea optimala a resurselor umane). Fie un numarde n persoane disponibile pentru m sarcini. Valoarea atribuita unei zilede lucru a persoanei i la o sarcina j este definita de cij , pentru i = 1, . . . , nsi j = 1, . . . , m. Problema se reduce la determinarea repartizarii optimede sarcini ce maximizeaza valoarea lucrului total. Repartizarea sarcinilorconsta n alegerea valorilor xij , ce reprezinta proportia din timpulpersoanei i consumata pe sarcina j pentru i = 1, . . . , n si j = 1, . . . , m.Vom avea astfel urmatoarele constrangeri:

mj=1

xij 1,ni=1

xij 1, xij 0.

Prima relatie denota ca o persoana nu poate depasi 100% din timpul saulucrand. A doua relatie precizeaza ca unei singure persoane i este permissa lucreze la o sarcina la un moment dat. In final, a treia relatie interzice


ca timpul petrecut pe o sarcina sa fie negativ. In concluzie, problemaenuntata anterior are urmatoarea formulare:

minxRnm

f(x)

(=

ni=1

mj=1

cijxij

)

s.l.: Cx d

mj=1

xij 1 i = 1, . . . , nni=1

xij 1 j = 1, . . . , mxij 0 i, j

.

Fie matricea coeficientilor cij data de

1 1 00.1 5 20 1 0.5

, adica n=m=3.

Sa se rezolve problema de mai nainte cu ajutorul functiei linprog.

Rezolvare. Initial compunem matricea C si vectorul d ce formeazaconstrangerile pentru a putea apela functia linprog n forma standard:

function [x,f]=alocare()

c=[1 -1 0 0.1 -5 2 0 -1 -0.5];

A1=[ones(1,3) zeros(1,6);...

zeros(1,3) ones(1,3) zeros(1,3);...

zeros(1,6) ones(1,3)];

A2=[eye(3) eye(3) eye(3)];

C = [A1; A2; -eye(9)]; d= [ones(6,1);zeros(9,1)];

[x f] = linprog(c,C,d);

end

Exemplul 9 (Proiectare filtre). Una dintre cele mai importanteprobleme ale ingineriei o reprezinta proiectarea filtrelor cu raspuns laimpuls prestabilit. Prin notiunea de filtru ntelegem un sistem dinamicliniar invariant n timp, care de cele mai multe ori este definit de relatiade convolutie dintre doua semnale discrete/continue:

y(t) =

u(s)h(t s)ds y(k) =

i=u(i)h(k i),


n care h este raspunsul la impuls corespunzator sistemului. Echivalent,n domeniul frecvential, transferul intrare-iesire pentru cazul continuuanterior se reprezinta sub forma:

Y () = H()U(), <


liniare de inegalitate, i.e. un program liniar:

minxRn,t

f(x, t) (= t)

s.l.: C

[xt

] d

(|T (i)

n1j=0

xjHj(i)| t i = 1, . . . , m).

Exercitiul 1. Se considera cazul particular ce cuprinde functiile: T () =cos2 + 3, H0() =

2

1+, H1() =

21+

, H2() = 2 + + 1 si vectorul

= [0.1 0.5 1 5 10]. Tinand cont ca valorile i reprezinta componentelevectorului , sa se scrie modelul problemei anterioare sub forma standarda unei probleme de programare liniara si sa se rezolve problema cuajutorul functiei Matlab linprog.

1.2.4 Programare patratica

Functia quadprog rezolva probleme de programare patratica, ce vizeazaminimizarea unei functii obiectiv patratice, n prezenta unui set deconstrangeri liniare. O problema de programare patratica se formuleazan general dupa cum urmeaza:

minxRn

1

2xTQx+ qTx

s.l.: Cx d, l x u,Ax = b,

unde Q = QT Rnn, q, l, u Rn, d Rm, b Rp, C Rmn, A Rpn.Sintaxa functiei quadprog cuprinde:

x = quadprog(Q,q,C,d,A,b,l,u,x0,optiuni)

[x,fval,exitflag,output,lambda] = quadprog(. . . ).

Pentru mai multe detalii, se poate apela comanda: help quadprog

Exemplul 10 (Problema Google). Ca urmare a progreselor tehnologicerecente, motoarele de cautare (e.g. Google, Yahoo) au devenit unpunct central de interes n optimizare, urmarind dezvoltarea de algoritmieficienti ce executa o cautare cat mai rapida. Constatand ca o cautareeficienta presupune o clasificare cat mai strategica a paginilor web,rezulta o reducere a problemei de cautare la una de clasificare a acestor


pagini. Problema de clasificare (ranking) se formuleaza usor n termenide grafuri ponderate (Fig. 1.7). De aceea, paginile web sunt interpretateca noduri ale unui graf, muchiile ca legaturi aferente dintre ele, iar fiecareimuchii i este asociata o pondere.

Figura 1.7: Exemplu graf orientat

Daca notam cu E matricea de adiacenta a grafului, elementele acesteiarespecta urmatoarea regula: Eij este nenul daca exista o muchiedirectionata de la i la j, altfel este nul. Daca tinem cont de faptul cagraful este ponderat si ca suma ponderilor este egala cu 1, atunci matriceaE va avea suma pe linii egala cu 1. Matematic, problema precedentapresupune gasirea unui vector propriu x corespunzator valorii proprii 1,i.e. sa se determine x Rn astfel ncat Ex = x. Reformuland maideparte problema n termeni de optimizare avem:

minxRn

f(x)

(=

1

2Ex x2

)s.l.: Cx d (x 0)

Ax = b

(ni=1

xi = 1

).

De cele mai multe ori, matricea are dimensiuni foarte mari (e.g. 109109)ceea ce face problema foarte greu de rezolvat (vezi Sectiunea 7.2.3).

Exercitiul 2. Fie matricea E =

0.1 0 0 0.90 0 0.5 0.50.2 0 0.2 0.60 0 0 1

. Sa se rezolve

problema cu ajutorul functiei quadprog, folosind mai ntai doarconstrangeri de egalitate si apoi incluzand si constrangerile de inegalitate.Ce se observa?

1.3. Probleme rezolvate de seminar 29

1.3 Probleme rezolvate de seminar

1.3.1 Diferentiere multivariabila

Problema 1. Fie functiile f : R3 R si g : R2 R definite de:

f(x) = x31x2 + 2x23x1 x2x3,

g(x) = ex1x2 + e2x11 e2x22.

(i) Sa se determine expresiile gradientilor f(x),g(x).(ii) Sa se determine expresiile Hessianelor 2f(x),2g(x).

Rezolvare. Pentru o functie diferentiabila de 2 ori f : Rn R, reamintimdefinitiile gradientului si Hessianei:

f(x) =

f(x)x1...

f(x)xn

, 2f(x) =

2f(x)

x21

. . . 2f(x)

x1xn...

...2f(x)xnx1

. . . 2f(x)x2n

.

(i) Expresiile gradientilor functiilor din enunt au forma:

f(x) =3x21x2 + 2x23x31 x34x3x1 x2

, g(x) = [ ex1x2 + 2e2x11ex1x2 2e2x22

].

(ii) Expresiile Hessianelor functiilor din enunt au forma:

2f(x)=6x1x2 3x21 4x33x21 0 1

4x3 1 4x1

,

2g(x)=[ex1x2+4e2x11 ex1x2ex1x2 ex1x24e2x22

].

Problema 2. Fie functia multidimensionala h : R3 R2,

h(x) =

[x22 + x3x21 + x2

].

(i) Determinati expresia Jacobianului h(x).


(ii) Aratati ca, n orice vector x R3, Jacobianul calculat are liniileliniar independente.

Rezolvare. (i) Din definitia Jacobianului pentru o functie h : R3 R2,i.e. h(x) =

[h1(x)x1

h1(x)x2

h1(x)x3

h2(x)x1

h2(x)x2

h2(x)x3

.

], avem h(x) =

[0 2x2 12x1 1 0

].

(ii) Din definitia proprietatii de liniar independenta avem: doi vectoriu, v Rn sunt liniar independenti daca nu exista un scalar 6= 0 astfelncat u = v. Aceasta proprietate este evidenta observand ca egalitatea:

2x110

=

02x2

1

nu poate avea loc daca 6= 0.

1.3.2 Probleme de optimizare generale

Problema 3. Fie urmatoarea problema de optimizare:

minxR2

x31 + x32

s.l.: (x1 + x2)2 0,

x1x32 = 1.

(i) Sa se precizeze care sunt constrangerile de egalitate si care suntcele de inegalitate.

(ii) Sa se determine analitic punctul de optim si valoarea optima.

Rezolvare. (i) Se observa ca inegalitatea (x1 + x2)2 0 este echivalenta

cu (x1+x2)2 = 0. Concluzionam ca nu exista constrangeri de inegalitate

si avem doua constrangeri de egalitate: h1(x) = x1 + x2 = 0 si h2(x) =x1x

32 + 1 = 0.

(ii) Din punctul (i) rezulta ca x1 + x2 = 0 si facand sistem cu cea de-adoua constrangere, concluzionam ca multimea fezabila este formata dindoua puncte, anume (x1, x2) = (1,1) si (x1, x2) = (1, 1) . Punctele deoptim global vor fi (x1, x

2) = (1,1) si (x1, x2) = (1, 1), din moment

ce valoarea optima f = 0 este aceeasi pentru ambele.

Problema 4. Fie functia f : n x, f(x) = ln xTAx, unde n esteo multime, numita multimea simplex, si este data de n = {x Rn :


ni=1

xi = 1, xi 0 i = 1, . . . , n}, A Rnn este o matrice simetrica sipozitiva (i.e. toate elementele matricei sunt ne-negative).

(i) Sa se determine expresia gradientului f(x) si a Hessianei 2f(x).(ii) Sa se determine constanta Lipschitz a gradientului functiei f .

Rezolvare: (i) Prin calcule simple, obtinem ca expresia gradientului este

data de f(x) = 2AxxTAx

, iar a hessianei 2f(x) = 2AxTAx(2Ax)(2Ax)T(xTAx)2

=2A

xTAx 4Ax

xTAx

(Ax

xTAx

)T.

(ii) Gradientul este continuu n sens Lipschitz daca urmatoareainegalitate are loc:

f(x)f(y)2 Lx y2 x, y n.Echivalent, aceasta proprietate se exprima si n termenii matriceiHessiane:

2f(x)2 L x n.Obervam ca matricea (2Ax)(2Ax)

T

(xTAx)2este pozitiv semidefinita, de aceea avem

urmatoarea marginire a normei Hessianei:

2f(x) 2AxTAx

x n.Mai mult, observam ca

minxn

xTAx = minxn

i,j

Aijxixj minxn

(miniAii

)x2 = 1

N

(miniAii

).

De aici obtinem o aproximare a constantei Lipschitz corespunzatoaregradientului functiei din enunt:

2f(x) N 2AminiAii = L.

Problema 5. Fie P Rnn o matrice simetrica si inversabila. Sa searate ca functia patratica:

f(x) =1

2xTPx+ xT b

are valoare minima (i.e. este marginita inferior) daca si numai dacaP 0 (i.e. este pozitiv semidefinita). Mai mult, n cazul n care estemarginita, aratati ca punctul de minim x este unic si determinat deexpresia x = P1b.


Rezolvare. Observam ca

1

2

(x+ P1b

)TP(x+ P1b

)=

1

2xTPx+ bTx+

1

2bTP1b.

De aici rezulta,

f(x) =1

2xTPx+ bTx =

1

2

(x+ P1b

)TP(x+ P1b

) 12bTP1b.

Presupunem ca P are o valoare proprie negativa, e.g. (unde > 0)si notam vectorul propriu asociat valorii proprii negative cu u. Atunci,pentru orice R, 6= 0, considerand x = u P1b si tinand cont derelatia Pu = u, avem:

f(x) =1

2

(x+ P1b

)TP(x+ P1b

) 12bTP1b

=1

2uTPu 1

2bTP1b

= 122u22

1

2bTP1b.

In final, observam ca > 0 si avem libertatea de a alege oricat de mare,de unde rezulta ca functia nu are minimum (i.e. minimum se atinge la). Pentru ca functia sa aiba minimum este necesar ca P 0. Pentrua arata ultima parte a rezultatului, din reformularea lui f(x) anterioaraavem:

f(x) =1

2

(x+ P1b

)TP(x+ P1b

) 12bTP1b.

Observam ca 12(x+ P1b)T P (x+ P1b) 0, deci minimul functiei se

atinge atunci cand x+ P1b = 0.

Problema 6. Fie problema de optimizare :

minxRn

f(x)

(=

cTx+ d

eTx+ f

)s.l. Gx h, Ax = b,

cu domf(x) ={x Rn : eTx+ f > 0}. Aratati ca aceasta problema

poate fi adusa la forma de programare liniara.


Rezolvare. Notam u(x) = eTx + f > 0. Impartind relatia prin u(x),

aceasta devine: eT xu(x)

+ f 1u(x)

= 1 > 0. De asemenea, folosim notatia

y(x) = xu(x)

Rn si z(x) = 1u(x)

R. Reformuland functia obiectiv avem:

f(x) =cTx+ d

u(x)= cT

x

u(x)+ d

1

u(x)= cTy(x) + dz(x).

Folosind schimbarea de variabila precedenta, putem aduce problema deminimizare dupa x la una ce presupune minimizarea dupa y(x) si z(x).Pentru a schimba variabila constrangerilor, mpartim prin u ultimeledoua seturi de egalitati/inegalitati. Astfel, obtinem un LP:

minyRn,zR

cTy + dz

s.l.: eTy + fz = 1, Ay = bz

Gy hz.

Problema 7. Fie urmatoarea problema de optimizare:

minxR2

2x21 + x22 x1x2

s.l.: 2xTx1 2, (aTx b)2 0

5|x1| 1, 3x1+x2 = 1,

unde a R2 si b R sunt dati. Sa se arate ca toate functiile care descriuproblema de optimizare sunt fie liniare, fie patratice. Sa se precizeze caresunt constrangerile de egalitate si care sunt cele de inegalitate.Rezolvare. Se observa ca functia obiectiv f(x) = 2x21 + x

22 x1x2 este

patratica, i.e.

f(x) =1

2

[x1x2

]T [4 11 2

]T [x1x2

].

Constrangerea 2xTx1 2 este echivalenta cu xTx 1 1, i.e. xTx 2.

Deci prima constrangere este una de inegalitate descrisa de o functiepatratica g1(x) = x

Tx 2 0.Constrangerea (aTx b)2 0 este echivalenta cu aTx b = 0. Deci adoua constrangere este de egalitate descrisa de o functie liniara h1(x) =aTx b = 0.


Constrangerea 5|x1|1 1 este echivalenta cu |x1| 1 0. Deci a treiaconstrangere este de inegalitate, descrisa de doua functii liniare g2(x) =x1 0 0 si g3(x) = x1 1 0.Constrangerea 3x1+x2 = 1 este echivalenta cu x1 + x2 = 0. Deci a patraconstrangere este de egalitate, descrisa de o functie liniara h2(x) = x1 +x2 = 0.Problema 8. O functie se numeste monomiala daca se prezinta subforma:

f(x) = cxa11 xa22 ...x

ann ,

n care c > 0, ai R si x Rn++, i.e. xi > 0 pentru orice i = 1, . . . , n.Considerand problema de programare geometrica :

minxRn

f(x)

s.l.: gi(x) 1 i = 1, . . . , mhi(x) = 1 i = 1, . . . , p,

cu fi si hi functii monomiale. Aratati ca o astfel de problema deprogramare geometrica poate fi scrisa sub forma unui LP.

Rezolvare. Din teoria optimizarii stim ca punctul de optim al uneiprobleme de optimizare cu functia obiectiv f(x), este acelasi cu cel alproblemei cu functia obiectiv log f(x). Intr-adevar, n cazul neconstransconditiile suficiente de optimalitate pentru problema originala (cu functiaobiectiv f(x)) suntf(x) = 0. Observam ca pentru cazul compunerii cufunctia logaritm (cu functia obiectiv log f(x)) conditiile de optimalitate

se transforma n f(x)

f(x)= 0, sau echivalent n f(x) = 0. Precizam

ca pentru cazul constrans are loc o echivalenta logaritmica similara, i.eavand problema originala:

minxRn

cxa11 xa22 ...x

ann

s.l.: cixbi11 x

bi22 . . . x

binn 1 i = 1, . . . , m

ejxdj1

1 xdj2

2 . . . xdjnn = 1 j = 1, . . . , p,

compunand functia obiectiv si constrangerile cu functia logaritm,obtinem o problema de optimizare cu acelasi punct de optim:

minxRn

log(cxa11 xa22 ...x

ann )

s.l.: log(cixbi11 x

bi22 . . . x

binn ) 0 i = 1, . . . , m

log(ejxdj1

1 xdj2

2 . . . xdjnn ) = 0 j = 1, . . . , p.

1.4. Probleme propuse 35

Folosind schimbarea de variabila log xi = yi, observam urmatoareareformularea liniara a unei functii monomiale:

f(x) = log(cxa11 . . . xann ) = log c +

ni=1

ai log xi = log c+

ni=1

aiyi.

Astfel, rezulta o forma LP a problemei n variabila y.

Problema 9. Fie functia patratica f : R2 R, f(x) = x21x22+x1x2. Sase arate ca functia f are gradient Lipschitz n raport cu norma Euclidiana,i.e. exista L > 0 astfel ncat:

f(x)f(y)2 Lx y2, x, y R2.si sa se determine constanta Lipschitz corespunzatoare.

Rezolvare. Remarcam ca functia f se poate formula ca f(x) = 12xTQx,

unde Q =

[2 11 2

]. Verificand relatia de continuitate Lipschitz avem:

QxQy2 = Q(x y)2 Q2x y2, x, y R2.Inegalitatea reiese din definitia normei matriceale induse. Evident,constanta Lipschitz este data de L = Q2 = max(Q).

1.4 Probleme propuse

Problema 1. Sa se arate ca problema de optimizare (1.1) este problemapatratica (QP). Evidentiati matricea Hessiana si precizati daca este saunu pozitiv definita.

Problema 2. Fie functia neliniara g(x) = 5e2x + cosx 1. Sa seaproximeze cu un polinom de gradul 3 functia g dupa modelul descrisn Exemplul 3 pentru punctele y1 = 1, y2 = 1/2, y3 = 0, y4 = 1/2si y5 = 1. Sa se rezolve problema de optimizare cu ajutorul functiilorfminunc si quadprog.

Problema 3. Fie urmatoarea problema de optimizare constransa:

maxxRn

1

2Ax2

s.l.: x1 1,


unde A Rmn. Sa se reformuleze problema precedenta ca una patratica(QP) si sa se rezolve apoi aceasta formulare cu ajutorul functiei Matlabquadprog.

Problema 4. Sa se calculeze expresiile gradientului si Hessianeicorespunzatoare urmatoarelor functii:

(i) f(x) = x21 x22 2x23 + 2x2x3 + 3x1x2 (ii) f(x) = x31 + x2x23(iii) f(x) = ex1x2+ex2x3+ex3x1 (iv) f(x) = x1 ln x1+x2 ln x2+x3 ln x3.

Problema 5. Sa se arate ca functia patratica f(x) = x21 + 2x22 2x23

x2x3+3x1x2 are gradientul continuu n sens Lipschitz n raport cu normaEuclidiana. Sa se determine constanta Lipschitz L.

Problema 6. Se considera urmatoarea problema de optimizare:

minxRn

f(x1, x2)

s.l.: 2x1 + x2 1, x1 + 3x2 1x1 0, x2 0.

Sa se traseze grafic multimea fezabila. Pentru urmatoarele functiiobiectiv, sa se determine multimea optima si valoarea optima:(i) f(x1, x2) = x1 + x2. (ii) f(x1, x2) = x1 x2.(iii) f(x1, x2) = x1. (iv) f(x1, x2) = max{x1, x2}.(v) f(x1, x2) = x

21 + 9x

22.

Problema 7. Se considera problema de optimizare:

minxRn

(x1 4)2 + (x2 2)2

s.l.: 4x21 + 9x22 36, 2x1 3

x21 + 4x22 = 4.

(i) Sa se reprezinte grafic schema multimii fezabile, multimile izonivel alefunctiei obiectiv si identificati punctul de optim pe grafic.(ii) Sa se rezolve din nou punctul (i) unde n enuntul problemeiminimizarea se nlocuieste cu maximizarea.

Problema 8. O fabrica furnizeaza patru tipuri de produse. Unul dintrematerialele brute necesare pentru fabricatie se afla n cantitate redusa,fiind disponibila numai o cantitate R de material. Pretul de vanzare alprodusului i este Si per kilogram. Mai mult, fiecare kilogram de produs iconsuma o cantitate ai de material brut (cel n cantitate redusa). Costul


de productie a xi kilograme din produsul i, excluzand costul materialuluin cantitate redusa, este de kix

2i , unde ki > 0 este cunoscut. Sa se dezvolte

un model matematic de optimizare pentru problema si apoi sa se rezolvecu una din functiile Matlab precizate anterior.

Problema 9. Se presupune ca cererile d1, . . . , dn pentru un anumitprodus n intervalul a n perioade de timp sunt cunoscute. Cererea ntimpul perioadei j poate fi satisfacuta din cantitatea produsa xj peparcursul perioadei sau din stocul magaziei. Orice exces de productiepoate fi stocat n magazie. Cu toate acestea, magazia are o capacitatelimitata K, iar costul stocarii unei unitati de la o perioada la alta este c.Costul productiei pe parcursul perioadei j este f(xj), pentru j = 1, . . . , n.Daca stocul initial este I0, sa se formuleze problema ca un programneliniar.

Problema 10. Intr-o localitate se urmareste amplasarea optima a unuinumar de depozite n vecinatatea magazinelor. Fie un numar de nmagazine, cu pozitii si cereri cunoscute. Magazinul i se afla n pozitia(ai, bi) si dispune de o cerere (grad de solicitare a produselor) ri. Cererilevor fi satisfacute cu ajutorul a m depozite cu capacitati cunoscute.Folosim urmatoarele notatii: (xi, yi) reprezinta pozitia necunoscuta adepozitului si ci este capacitatea depozitului i pentru orice i = 1, . . . , m;dij este distanta de la depozitul i la magazinul j si wij reprezinta unitatiletransportate de la depozitul i la magazinul j pentru orice i = 1, . . . , m sij = 1, . . . , n.Sa se modeleze problema de gasire a pozitiilor optime corespunzatoaredepozitelor. Pentru masurarea distantelor n plan, se poate folosi normaEuclidiana (norma 2) sau orice alta norma, cu conditia sa fie precizata.

Problema 11. Fie matricele X, Y, Z Rnn. Sa se arate egalitatea:

Tr(XY Z) = Tr(ZXY ) = Tr(Y ZX).

Problema 12. Fie matricea simetrica Q Rnn. Sa se determineexpresia valorilor proprii extreme min si max sub forma valorilor optimecorespunzatoare unor probleme de optimizare.

Problema 13. Fie x Rn, pentru norma , definim norma duala:x = min

y1x, y. Sa se demonstreze urmatoarele relatii, pentru orice

x Rn:x1 = x, x = x1, x2 = x2.

Capitolul 2

Probleme de optimizare

convexa

2.1 Preliminarii

Fie problema de optimizare,

minxRn

f(x)

s.l.: gi(x) 0 i = 1, . . . , m (2.1)Ax = b.

Daca functia obiectiv f si functiile ce definesc constrangerile deinegalitate gi sunt convexe atunci problema se numeste problema deoptimizare convexa. Convexitatea detine un rol crucial n optimizare,deoarece problemele cu aceasta proprietate prezinta trasaturi teoreticefoarte bune (e.g., punctele de optim locale sunt, de asemenea, puncte deoptim globale); si ceea ce este mai important, pot fi rezolvate numericn mod eficient, ceea ce nu este valabil pentru problemele neconvexe.Vom studia mai departe notiunile de functie convexa si multime convexapentru a sublinia proprietatile remarcabile ale problemelor convexe.

2.1.1 Multimi convexe

Pentru o expunere clara si concisa, introducem urmatoarele definitii:

Definitia 1. Multimea S Rn se numeste convexa daca pentru oricaredoua puncte x1, x2 S si un scalar [0, 1] avem x1 + (1 )x2 S,

2.1. Preliminarii 39

i.e. segmentul generat de oricare doua puncte din S este inclus n S(Fig. 2.1).

x1

x2

x1

x2

Figura 2.1: Exemplu de multime convexa (stanga) si multime neconvexa(dreapta).

Exemplul 11. Fie multimea Ln ={[

xt

] Rn+1 : x t

}, denumita

si conul Lorentz sau conul de nghetata. Sa se traseze graficul acesteimultimi n Matlab pentru n = 2.

Rezolvare. Figura conului (Fig. 2.2) poate fi trasata de urmatoareasecventa de cod:

function create_cone

[y1,y2]=meshgrid(-1:0.01:1,-1:0.01:1);

y3=sqrt(y1.^2+y2.^2); dom=[y3>1];

z3=y3; z3(dom)=inf;

figure(1);

hold on;

surf(y1,y2,z3); set(gca,FontSize,15);

hold off;

end

Figura 2.2: Con de ordinul II (con Lorentz) n R3 pentru n = 2.

40 Capitolul 2. Probleme de optimizare convexa

Studiul multimilor convexe si al proprietatilor acestora faciliteaza analizamultimilor fezabile aferente problemelor de tipul (2.1) si a sirurilorgenerate de algoritmi cu restrictii la aceste multimi. In continuare vomintroduce si analiza notiunea de functie convexa.

2.1.2 Functii convexe

Definitia 2. Functia f : Rn R se numeste convexa daca domeniul sauefectiv domf este o multime convexa si daca urmatoarea inegalitate:

f(x1 + (1 )x2) f(x1) + (1 )f(x2), (2.2)

este satisfacuta pentru orice x1, x2 domf si [0, 1].

1 0.5 0 0.5 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

f( x1+(1)x2)

f(x)=x2

f(x1)+(1) f(x2)

f(x2)

f(x1)

x1 x2

Figura 2.3: Exemplu de functie convexa f(x) = x2.

Inegalitatea (2.2) se poate generaliza la un numar p de puncte, purtandnumele de inegalitatea lui Jensen: f este o functie convexa daca si numaidaca:

f

(pi=1

ixi

)

pi=1

if(xi)

pentru orice xi domf si i [0, 1], cupi

i = 1.

2.2. Pachet CVX 41

Interpretarea geometrica a convexitatii este foarte simpla (Fig. 2.3).Pentru o functie convexa, fie doua puncte din domeniul sau x, y domf ,atunci valorile functiei evaluate n punctele din intervalul [x, y] sunt maimici sau egale decat cele aflate pe segmentul cu capetele (x, f(x)) si(y, f(y)).

Remarca 1. O functie este convexa daca si numai daca restrictiadomeniului sau la o dreapta (care intersecteaza domeniul) este deasemenea, convexa. Cu alte cuvinte, f este convexa daca si numai dacaoricare ar fi x domf si o directie d, functia g() = f(x + d) esteconvexa pe domeniul { R : x+d domf}. Aceasta proprietate estefoarte utila n problemele ce implica teste de convexitate ale functiilor.

Pentru a sublinia legatura dintre notiunile introduse anterior, facemurmatoarea observatie: orice multime M definita ca multimea subnivela unei anumite functii f , i.e. M = {x Rn : f(x) c}, unde c esteo constanta, este convexa daca si numai daca functia f este convexa(pentru demonstratie vezi [1]).

Exemplul 12. Fie functia convexa f : R2 R, definita de f(x) =12xT[2 11 0.5

]x+ [0.5 1]x. Sa se traseze graficul functiei n Matlab.

Rezolvare. Graficul functiei (Fig. 2.4) este generat de urmatoareasecventa de cod (unde am folosit functia gradient(z) pentru a descrieculoarea, graduala de-a lungul axei z, graficului rezultat):

function create_parab

[x,y] = meshgrid([-10:.2:10]);

z=x.^2+0.25*(y.^2)-x.*y-0.5*x+y;

surf(x,y,z,gradient(z))

end

2.2 Pachet CVX

In aceasta subsectiune expunem o scurta prezentare a pachetului deoptimizare CVX. CVX reprezinta un sistem de modelare a problemelorde optimizare convexa pe baza limbajului Matlab. CVX realizeazatransformarea comenzilor Matlab ntr-un limbaj de modelare, cu


105 0 5 10 10

50

510

50

0

50

100

150

200

250

Figura 2.4: Graficul functiei patraticef(x) = x21 + 1/4x

22 x1x2 1/2x1 + x2.

posibilitatea de declarare a constrangerilor si functiei obiectiv prinintermediul sintaxei Matlab standard.De exemplu, consideram urmatorul model de optimizare convexa:

minxRn

Ax b2s.l.: Cx = d, x 1.

Urmatoarea secventa de cod genereaza si rezolva o instanta aleatorie aacestui model:

m = 10; n = 8; p = 2;

A = randn(m,n); b = randn(m,1);

C = randn(p,n); d = randn(p,1);

cvx_begin

variable x(n)

minimize(norm(A*x - b, 2))

subject to

C*x == d

norm(x, Inf)

2.2. Pachet CVX 43

se porneste Matlab-ul; se comuta directorul curent n directorulunde am realizat dezarhivarea si se executa n consola comandacvx_setup.

Remarca 2. In unele cazuri (de cele mai multe ori n Linux) este nevoiede crearea sau modificarea fisierului startup.m ce elimina necesitateaca utilizatorul sa introduca comanda cvx_setup la fiecare pornire aMatlab-ului.

2.2.2 Elemente de baza

Daca instalarea este realizata cu succes, ncepem prin a preciza caorice program CVX se scrie n interiorul unei functii Matlab. Pentrua diferentia continutul codului CVX de restul programului Matlab,programele CVX se delimiteaza cu comenzile cvx_begin si cvx_end.Valorile variabilelor create n portiunea de cod Matlab se pot folosi caparametri n problemele de optimizare rezolvate cu CVX. La nceputuloricarui program CVX se definesc variabilele de decizie si dimensiunileacestora, e.g.:

m = 10; n = 8; A = randn(m,n); b = randn(m,1);

cvx_begin

variable x(n)

minimize( norm(A*x-b) )

cvx_end

In acest exemplu, variabila de decizie este vectorul x din Rn, iar matriceaA si vectorul b reprezinta variabile Matlab ce sunt utilizate ca parametrin cadrul CVX. Variabilele CVX se declara folosind comanda variable sispecificarea dimensiunii variabilei, nainte de a se utiliza n constrangerisau n expresia functiei obiectiv. Sintaxa declararii unei variabile poatefi una din urmatoarele forme:

variable x(10) variable Y(20,10) variable Z(5,5,5).

Sunt disponibile o varietate de optiuni aditionale pentru precizareastructurii matriceale, n cazul n care variabila de decizie este de tipmatrice cu proprietati speciale (simetrica, Toeplitz) e.g.:

variable Y(50,50) symmetric

variable Z(100,100) hermitian toeplitz.


Pentru lista ntreaga a optiunilor se poate consulta adresa:

http://cvxr.com/cvx/doc/basics.html.

Declararea functiei obiectiv a problemei de optimizare necesita precizareatipului de problema (e.g. minimizare, maximizare) prin intermediulcuvintelor cheie minimize si maximize, e.g.

minimize( norm( x, 1 ) )

maximize( geo_mean( x ) ).

Precizam ca este imperativ ca functia obiectiv sa fie convexa cand folosimminimize si concava cand folosim maximize. In caz contrar, pachetulCVX va furniza un mesaj de eroare corespunzator. Constrangerilesuportate de modelele CVX sunt cele de egalitate (liniare) impuse prinoperatorul == si de inegalitate impuse de operatorii =. Pentruconstrangeri de tip box (lant de inegalitati) este disponibila sintaxal


1.5 1 0.5 0 0.51

0.8

0.6

0.4

0.2

0

0.2

0.4

0.6

0.8

1

x1x 2

Bila de raza maxima inscrisa intrun poliedru 2D

Figura 2.5: Bila de raza maxima nscrisa ntr-un poliedru.

variable x_c(2)

maximize ( r )

a1*x_c + r*norm(a1,2)


Rezolvare. Formularea problemei precedente n termeni de optimizareconduce la rezolvarea urmatoarei probleme:

minH

max|H()Hd()|, (2.3)

unde H este functia raspunsului n frecventa, iar variabila de decizieare rolul raspunsului la impuls. O aproximare convenabila a problemeianterioare presupune extragerea unui set finit m de frecvente i, cu i =1, . . . , m, si definirea unei matrice A:

A =

1 ej1 e2j1 . . . enj1

1 ej2 e2j2 . . . enj2

. . . . . . . . . . . .1 ejm e2jm . . . enjm.

rezultand o aproximare a problemei (2.3), i.e.

minxRn

max1im

|AixHd(i)|, (2.4)

unde Ai reprezinta linia i a matricei A.Mai departe, prezentam secventa de cod ce rezolva problema deoptimizare (2.4) si determinam aproximarea optima a functiei de transferHd() = e

jD (reprezentata n Fig. 2.6):

0 5 10 15 200.2

0

0.2

0.4

0.6

0.8

1

1.2

n

h(n)

Figura 2.6: Aproximarea optima n sens Chebyshev a functiei de transferHd() = e

jD.

n = 20; m = 15*n;

w = linspace(0,pi,m); % omega

% Construim un raspuns in frecventa dorit


D = 8.25; % valoare intarziere

Hdes = exp(-j*D*w); % raspuns in frecventa dorit

% Filtru Gaussian cu faza liniara (decomentati liniile

% de mai jos pentru proiectare)

% var = 0.05;

% Hdes = 1/(sqrt(2*pi*var))*exp(-(w-pi/2).^2/(2*var));

% Hdes = Hdes.*exp(-j*n/2*w);

% Rezolvam problema minimax de proiectare a filtrului

% A reprezinta matricea folosita la calculul

% raspunsului in frecventa

% A(w,:) = [1 exp(-j*w) exp(-j*2*w) ... exp(-j*n*w)]

A = exp( -j*kron(w,[0:n-1]) );

% formularea optimala a filtrului Chebyshev

cvx_begin

variable h(n,1)

minimize( max( abs( A*h - Hdes ) ) )

cvx_end

% verificam daca problema a fost rezolvata cu succes

disp([Problem is cvx_status])

if ~strfind(cvx_status,Solved)

h = [];

end

% Generam figurile aferente filtrului

figure(1); stem([0:n-1],h); xlabel(n); ylabel(h(n));

% figura raspunsului in frecventa

H = [exp(-j*kron(w,[0:n-1]))]*h;

figure(2)

% magnitudine

subplot(2,1,1);

plot(w,20*log10(abs(H)),w,20*log10(abs(Hdes)),--);

xlabel(w); ylabel(mag H in dB); axis([0 pi -30 10]);

legend(optimizat,dorit,Location,SouthEast)

% faza


subplot(2,1,2); plot(w,angle(H));

axis([0,pi,-pi,pi]); xlabel(w), ylabel(faza H(w));

2.4 Probleme rezolvate de seminar

2.4.1 Multimi si functii convexe

Problema 1. Fie multimea S descrisa de:

S = {x Rn : aTi x = bi, cTj x dj, i = 1, . . . , m, j = 1, . . . , p}= {x Rn : Ax = b, Cx d}.

Sa se demonstreze ca multimea este convexa (multimile definite deegalitati si inegalitati liniare, cum este cea din enunt, se numesc poliedre).

Rezolvare. Se arata ca pentru orice x1, x2 S si [0, 1] rezulta x1 +(1 )x2 S. Daca x1, x2 S atunci avem:{

Ax1 = b, Cx1 dAx2 = b, Cx2 d.

Notand x = x1 + (1 )x2, putem finaliza demonstratia prinurmatoarea observatie:

Ax = Ax1 + (1 )Ax2 = b+ (1 )b = bCx = Cx1 + (1 )Cx2 d+ (1 )d = d.

In acest fel am demonstrat ca multimea definita anterior este convexa.

Problema 2. Sa se demonstreze ca urmatoarele multimi se pot definisub forma unor poliedre:

(i) S1 = {x Rn : x 0, x1 = 1}

(ii) S2 = {x Rn : x 1}

(iii) S3 = {x Rn : x1 1}.


Rezolvare. : (i) In cazul multimii S1 observam urmatoarele echivalente:

S1 =

{x Rn : x 0,

ni=1

xi = 1

}=

{x Rn : x 0,

ni=1

xi = 1

}

= {x Rn : Inx 0, [1 . . . 1]x = 1} ,unde In reprezinta matricea identitate de ordin n. Ultima formularedenota forma poliedrala a multimii S1.(ii) De asemenea, n cazul multimii S2, urmand acelasi rationament:

S2 =

{x Rn : max

1in|xi| 1

}= {x Rn : |xi| 1, i = 1, . . . , n}

= {x Rn : 1 xi 1 i = 1, . . . , n}= {x Rn : xi 1, xi 1, i = 1, . . . , n}

=

x Rn :

[InIn

]x

1...1

,

obtinem o forma poliedrala.(iii) Pentru a determina forma poliedrica a multimii S3 definim o multimepoliedrica auxiliara:

Q =

{[xT tT ]T R2n :

ni=1

ti = 1, |xi| ti}.

Mai departe, definim proiectia unei multimi poliedrice pe un subspatiu dedimensiuni reduse si folosim aceasta notiune pentru a studia multimea S3.

Definitia 3. Fie multimea P = {x Rn1, y Rn2 : Ax + By b},unde A Rmn1 , B Rmn2 , b Rm. Proiectia multimii P pe subspatiulvariabilelor x este data de:

Px = {x Rn1 : y Rn2 , [xT yT ]T P}.Proiectia multimii Q pe subspatiul variabilelor x conduce la urmatoareaserie de echivalente:

Qx = {x Rn : t Rn a..ni=1

ti = 1, |xi| ti}

= {x Rn :ni=1

|xi| 1} = {x Rn : x1 1} = S3.


Tinand cont ca proiectia oricarei multimi poliedrice este o multime polie-drica, concluzionam ca multimea S3 este de tip poliedric.

Problema 3. Fie multimile:

(i) Ln ={[

xt

] Rn+1 : x t

}numit si con Lorentz sau con de

ordinul II;

(ii) Sn+ = {X Sn : X 0} conul semidefinit.Sa se demonstreze ca multimile precedente sunt conuri auto-duale.

Rezolvare: (i) Din definitia conului dual avem:

Ln = {y Rn : x, y 0 x Ln}.

Problema se reduce la a demonstra ca Ln = Ln, ceea ce este echivalentcu satisfacerea simultana a incluziunilor: Ln Ln si Ln Ln.Aratam acum prima incluziune. Fie y =

[y1v

] Ln, atunci pentru

orice x =

[x1t

] Ln avem:

x, y = yT1 x1 + vt 0. (2.5)

Stiind ca x Ln atunci x t. Ramane sa demonstram ca y v.Din ipoteza ca (2.5) are loc pentru orice vector x Ln, atunci inegalitateaeste satisfacuta, de asemenea, pentru un x ales. Alegand x =

[x1t

]=[ y1y1

1

]si nlocuind n (2.5) obtinem:

yT1 y1y1 + v = y1+ v 0,

din care rezulta prima incluziune.

Pentru a doua incluziune, fie y =

[y1v

] Ln, x =

[x1t

] Ln. Din

inegalitatea Cauchy-Schwartz |x, y| xy rezulta:

yT1 x1 + vt > y1x1+ vt vt + vt = 0,


unde n a doua inegalitate am utilizat ipoteza ca x, y Ln. In concluzie,pentru orice y Ln avem ca y Ln.(ii) In mod similar, aratam prin dubla incluziune ca Sn+ = Sn+ . Fie Y Sn+ , atunci Tr(Y X) 0. Din urmatoarele relatii:

xTY x = Tr(xTY x) = Tr(Y xxT ) = Tr(Y X) 0,

se deduce Y Sn+. Pentru a doua incluziune, presupunem Y,X matricepozitiv semidefinite. Remarcam urmatoarea relatie:

Tr(YX) = Tr(Y V TV ) = Tr(V Y V T) 0, (2.6)

n care am folosit descompunerea valorilor proprii corespunzatoarematricei X si proprietatea de permutare a functiei matriceale Tr().In concluzie, Y Sn+ datorita relatiei (2.6), de unde reiese a douaincluziune.

Problema 4. Sa se demonstreze ca urmatoarele functii sunt convexe pemultimile specificate:(i) f(x) = log x, domf = (0,).(ii) f(x) = 1

2xTQx+ qTx+ r, domf = Rn, Q 0.

(iii) f(x, t) = xTxt, domf = Rn (0,).

(iv) f(x) = Ax b, A Rmn, domf = Rn.(v) f(X) = log detX, domf = Sn++.Rezolvare: (i) Observam ca functia f(x) = log x satisface conditiile deordinul II ale convexitatii:

f (x) =1

x2> 0 x > 0.

(ii) Hessiana functiei f(x) = 12xTQx+ qTx+ r este data de 2f(x) = Q.

Observam ca functia este convexa deoarece matricea Q 0 (este pozitivsemidefinita).

(iii) In aceeasi maniera aratam ca functia f(x, t) = xTxt

este convexa pedomeniul Rn (0,]. Din definitia Hessianei rezulta:

2f(x) =[

2tIn 2xt

2xTt2

2xTxt3

],

unde cu In am notat matricea identitate de ordin n. Pentru a determinadaca functia ndeplineste conditiile de ordinul II ale convexitatii observam


ca:

[uTvT ]2f(x)[uu

]= [uTvT ]

[2tu 2xv

t2

2xT ut2

+ 2xTxvt3

]

=2

tuTu 2u

Txv

t2 2x

Tuv

t2+

2xTxv2

t3

=2

t3(t2uTu 2tuTxv + xTxv2)

=2

t3tu xv2.

Tinand cont de domeniul de definitie al functiei, remarcam ca pentru

orice vector

[uv

], termenul drept din ultima egalitate este pozitiv. De

aici este evidenta proprietatea de pozitiv definire a matricei Hessianecorespunzatoare functiei.(iv) Deoarece functia este nediferentiabila n punctul 0, observam cafunctia f(x) = Ax b nu este diferentiabila n punctele x ce satisfacAx = b. Fie doua puncte din multimea domf pentru care f(x1) =Ax1 b, f(x2) = Ax2 b. Pentru a arata convexitatea functiei f ,notam x = x1 + (1 )x2 si deducem sirul de relatii:

f(x) = f(x1 + (1 )x2)= A (x1 + (1 )x2) b (Ax1 b)+ (1 )(Ax2 b)= Ax1 b + (1 )Ax2 b= f(x1) + (1 )f(x2).

(v) Fie functia f(X) = log detX, X Sn++. Aratam convexitatea luif prin intermediul reducerii domeniului acesteia la o dreapta. Mai exact,folosim urmatoarea proprietate a functiilor convexe: f este convexadaca functia scalara ce se obtine din restrictionarea la o dreapta este deasemenea convexa (conform Remarca 1). Revenind la functia matriceala,consideram matricele X Sn++, D Sn si aratam ca functia scalarag(t) = f(X + tD) este convexa, observand urmatoarele egalitati:

g(t) = log det(X + tD)= log det(X1/2X1/2 + tD)= log det (X1/2 (In + tX1/2DX1/2)X1/2)


= log (detX1/2 det (In + tX1/2DX1/2) detX1/2)= log (detX det (In + tX1/2DX1/2))= log detX log det (In + tX1/2DX1/2) .

Pe de alta parte, stiind ca pentru orice matrice A Sn, cu spectrul(A) = {1, . . . , n}, transformarea B = In + tA, t R, modificaspectrul astfel ncat (B) = {1+ t1, . . . , 1+ tn} si detB =

ni=1

(1+ ti).

Pentru a aplica aceasta proprietate n sirul de relatii precedent, notamZ = X1/2DX1/2 avand spectrul (Z) = {1, . . . , n}. Rescriind g(t),rezulta:

g(t) = log detX logni=1

(1 + ti)

= log detX ni=1

log(1 + ti).

Mai departe, observam ca functia h(t) = log(1 + tu) din componentacelei anterioare, este convexa, deoarece cea de-a doua derivata satisface:

h(t) =u

(1 + tu)2 0 u R.

Stiind ca functia definita de o suma de functii convexe este convexa,ajungem la concluzia ca g(t) este convexa n t. Deci, f(X) este convexa.

Problema 5. Fie functia f : Rn R. Sa se determine functia conjugataf (y) pentru urmatoarele exemple:

(i) f(x) = ex (ii) f(x) = x log x

(iii)f(x) = 12xTQx,Q 0 (iv) f(x) = log

ni=1

exi .

Rezolvare. Definim functia conjugata corespunzatoare functiei f :

f (y) = maxxdomf

y, x f(x).

In primele doua cazuri monovariabile (i) si (ii), observam ca functiaobiectiv a problemei de maximizare este concava pe domeniul functiei


f , de aceea conditiile de optimalitate de ordinul I conduc la urmatoareleexpresii:

(i) f (y) = y log y y, (ii) f (y) = ey1.(iii) In acest caz, conjugata functiei f are forma f (y) = max

xRnyTx

12xTQx. Datorita proprietatii de convexitate, observam ca solutia

problemei de optimizare este data de x = Q1y. Inlocuind n expresiafunctiei conjugate rezulta f (y) = yTQ1y 1

2yTQ1y = 1

2yTQ1y.

(iv) Conditiile necesare de optimalitate de ordinul I corespunzatoare

problemei de maximizare maxxRn

y, x logni=1

exi se reduc la urmatoarele

relatii:

ex

i

ni=1

ex

i

= yi, i = 1, . . . , n. (2.7)

Observam ca relatia (2.7) este satisfacuta si functia f (y) ia valori finite

numai n cazul n care argumentul functiei conjugate satisfaceni=1

yi =

1, y 0. Presupunand ca argumentul y satisface cele doua conditii, avemxi = log

ni=1

ex

i + log yi pentru orice i = 1, . . . , n. Substituind n expresia

functiei conjugate rezulta:

y, x logni=1

ex

i=

ni=1

yi

(log

ni=1

ex

i + log yi

) log

ni=1

ex

i

=ni=1

yi log yi + logni=1

ex

i

(ni=1

yi 1)=

ni=1

yi log yi.

In concluzie, functia conjugata f (y) este definita de:

f (y) =

ni=1

yi log yi, dacani=1

yi = 1, y 0,, altfel.


Problema 6. Sa se demonstreze ca urmatoarea problema de optimizareeste convexa:

minxR2

f(x) = x12 + x2

2

s.l.: g1(x) =x1

1 + x22 0, g2(x) = ex1+x2 1 0

h(x) = (x1 x2 1)2 = 0.Rezolvare. Pentru a arata convexitatea problemei din enunt:

minxR2

f(x)

s.l.: g1(x) 0, g2(x) 0, h(x) = 0.este suficient sa demonstram convexitatea functiilor f , g1 si g2, siliniaritatea functiei h. Observam ca Hessiana functiei f are formaexplicita 2f(x) = I2 0, unde I2 este matricea identitate de ordinul II.Deci, functia f este convexa deoarece satisface conditia de convexitate deordinul II. In cazul functiei g1, distingem faptul ca inegalitatea

x11+x22

0este satisfacuta doar n cazul n care x1 0. In concluzie, constrangereag1(x) 0 este echivalenta cu o constrangere liniara (si deci convexa)x1 0. Pentru inegalitate g2(x) 0, observam ca este echivalentacu x1 + x2 0, i.e. este si aceasta echivalenta cu o constrangereliniara. Similar, pentru egalitatea definita de functia h, gasim urmatoareaechivalenta ntre multimi:

{x R2 : h(x) = 0} = {x R2 : x1 x2 = 1},rezultand o functie liniara n x. In final, putem rescrie problema subforma unui QP convex:

minxR2

x12 + x2

2

s.l.: x1 0, x1 + x2 0, x1 x2 = 1.

Problema 7. Sa se demonstreze ca urmatoarea problema de optimizareeste convexa:

minxR2

f(x) (= log(aTx b))

s.l.: g1(x) = exT x e 0, g2(x) = (cTx d)2 1 0

h(x) = (x1 + 2x2)4 = 0.


Rezolvare. Pentru a demonstra ca o problema de optimizare este convexa,trebuie sa demonstram ca functia obiectiv este convexa si multimeafezabila definita de constrangeri este convexa. Pentru functia f(x) = log(aTx b) deducem expresia gradientului si a Hessianei:

f(x) = 1aTx ba,

2f(x) =1

(aTx b)2aaT

unde avem desigur ca (aTx b)2 > 0. Daca notam y = aTx, observamca:

xTaaTx = yTy = y22 0,deci matricea aaT este pozitiv semidefinita. Drept rezultat, 2f(x) 0peste ntregul dom f = {x Rn : aTx b > 0}, i.e. satisfaceconditiile de convexitate de ordin II si deci functia obiectiv este convexa.Pentru a arata ca multimea constrangerilor este convexa, este suficient saaratam ca functiile din constrangerile de egalitate sunt liniare, iar cele dinconstrangerile de inegalitate sunt convexe. Observam ca constrangereah(x) = 0 este echivalenta cu egalitatea x1 + 2x2 = 0, a carei functie esteliniara. Pentru g1(x) 0, observam ca este echivalenta cu xTx 1 0.Functia xTx 1 este o functie patratica, diferentiabila de doua ori, cuHessiana 2I2 0, deci satisface conditiile de convexitate de ordin II sieste implicit convexa. Constrangerea g2(x) 0 este echivalenta cu:

1 (cTx d) 1{cTx d 1cTx+ d 1

[cT

cT]x+

[dd

][11

],

deci aceasta se reduce la o constrangere de inegalitate unde functia esteliniara si implicit convexa.

Problema 8. Sa se determine problema convexa de programare semidefi-nita ce aproximeaza urmatoarea problema neconvexa:

maxxRn

xTATAx

s.l.: x2 1,

unde A Rmn si x2 =

ni=1

x2i .


Rezolvare. Reamintim ca pentru orice Q Rnn si x Rn, functiaTr(Q) satisface relatia: Tr(xTQx) = Tr(QxxT ). Pe baza acestei relatii,problema precedenta se scrie sub urmatoarea forma echivalenta:

maxXRnn

Tr(ATAX

)s.l.: rang(X) = 1, T r(X) = 1.

Se obtine relaxarea convexa prin renuntarea la constrangerea de egalitateneliniara rang(X) = 1. In concluzie, avem urmatoarea aproximareconvexa a problemei de optimizare originala:

maxXRnn

Tr(ATAX

)s.l.: Tr(X) = 1.

2.5 Probleme propuse

Problema 1. Sa se determine care dintre urmatoarele functii suntconvexe, concave sau niciuna dintre cele doua variante. Sa seargumenteze rezultatele obtinute.

(i) f(x1, x2) = 2x21 4x1x2 8x1 + 3x2,

(ii) f(x1, x2) = x1e(x1+3x2),

(iii) f(x1, x2) = x21 3x22 + 4x1x2 + 10x1 10x2,(iv) f(x1, x2, x3) = 2x1x2 + 2x

21 + x

22 + 2x

23 5x1x3,

(v) f(x1, x2, x3) = 2x21 3x22 2x23 + 8x1x2 + 3x1x3 + 4x2x3.Problema 2. Sa se determine submultimea din {x R : x > 0} pe carefunctia f(x) = eax

b

este convexa. Parametrii a, b satisfac a > 0, b 1.

Problema 3. Sa se determine n domeniul axei reale n care functiaf(x) = x2(x2 1) este convexa.

Problema 4. Fie functiile convexe f1, . . . , fm : Rn R. Sa se

demonstreze ca urmatoarele compuneri ale acestora sunt convexe:


(i) g(x) =mi=1

ifi(x), i > 0,

(ii) h(x) = max{f1(x), . . . , fm(x)}.

Problema 5. Fie multimea K = {x Rn : x 0}. Sa se demonstrezeca multimea K este un con si sa se determine conul dual corespunzator.

Problema 6. Sa se demonstreze ca functia f(x) =ri=1

ix[i] este convexa

n x, unde 1 r 0, iar x[n] reprezinta a n-a cea mai marecomponenta din vectorul x.

Problema 7. Determinati functiile conjugate corespunzatoare functiilor:

(i) f(x) = max1in

xi cu domeniul Rn;

(ii) f(x) = xp cu domeniul R++, unde p > 1. Comentati cazul n carep < 0.

(iii) f(x) =

(

ni=1

xi

)1/ncu domeniul Rn++;

(iv) f(x, t) = log(t2xTx) cu domeniul {(x, t) Rn R : x2 t}.Problema 8. Sa se demonstreze ca functia:

f(x) =1

x1 1

x2 1

x3 1

x4

,

definita pe domeniul n care toti numitorii sunt pozitivi, este convexa sistrict descrescatoare. (Cazul cand n = 4 nu prezinta nicio particularitate,caracteristicile functiei se mentin si pentru cazul general n oarecare.)

Problema 9. Presupunem ca functia f este convexa, iar 1 > 0, i 0pentru i = 2, . . . , n si

ni=1

i = 1. Fie punctele x1, . . . , xn domf . Sa sedemonstreze ca are loc inegalitatea:

f(1x1 + + nxn) 1f(x1) + + nf(xn).


Problema 10. Sa se arate ca o functie f : R R este convexa daca sinumai daca domeniul domf este convex si

det

1 1 1x y zf(x) f(y) f(z)

0,

pentru orice x, y, z domf si x < y < z.Problema 11. Sa se arate ca maximul unei functii convexe pestepoliedrul P = conv {v1, . . . , vk} este atins ntr-unul dintre varfurilepoliedrului, i.e.,

supxP

f(x) = maxi=1,...,k

f(vi).

Indiciu: Se presupune ca afirmatia este falsa si se utilizeaza inegalitatealui Jensen.

Problema 12. Fie o functie convexa f si o functie g definita dupa cumurmeaza:

g(x) = min>0

f(x)

.

(i) Sa se demonstreze ca functia g este omogena, i.e. g(tx) = tg(x)pentru orice t 0.

(ii) Sa se demonstreze ca daca o functie h este omogena si h(x) f(x)pentru orice x, atunci avem h(x) g(x) pentru orice x.

(iii) Sa se demonstreze ca functia g definita anterior este convexa.

Problema 13. Fie functia omogena f : Rn R, i.e. f(tx) = tf(x)pentru orice t 0 si x Rn. O functie omogena se numeste subaditivadaca satisface, n plus, urmatoarea relatie:

f(x) + f(y) f(x+ y).Sa se demonstreze ca, pentru functiile omogene, subaditivitatea esteechivalenta cu convexitatea.

Problema 14. Fie o multime S nevida, marginita si convexa n Rn si ofunctie f definita de:

f(x) = maxyS

yTx.

Functia f se numeste functia suport a multimii S.


(i) Sa se arate ca functia f este omogena si convexa.

(i) Sa se determine explicit functia suport pentru multimile S definitede

S = {x Rn : Axp 1} p = 1, 2,,unde matricea A este inversabila.

Problema 15. Fie functia convexa f : Rn R. Definim proprietateade convexitate tare, enumerand conditiile de ordin 0, I si II (ce depindde gradul de diferentiabilitate al functiei), ntr-un mod similar cu celecorespunzatoare cazului convex. Pentru > 0, functia f se numestetare convexa n raport cu norma-p daca:

(Conditii de ordinul 0) functia f satisface urmatoarea inegalitatepentru orice x, y Rn:

f(x+ (1 )y) f(x) + (1 )f(y) (1 )2

x y2p.

(Conditii de ordinul I) functia f este diferentiabila si satisfaceurmatoarea inegalitate pentru orice x, y Rn:

f(y) f(x) + f(x), y x+ 2x y2p.

(Conditii de ordinul II) functia f este diferentiabila de doua ori sisatisface urmatoarea inegalitate pentru orice x Rn:

xT2f(x)x x2p,

unde xp =(

ni=1

xpi

) 1p

, 1 p .

Sa se demonstreze ca:

(i) Pentru norma 2, conditiile de ordin I implica conditiile de ordin 0,i.e. daca functia f satisface conditiile de ordin I atunci ea satisfacesi conditiile de ordin 0.

(ii) Pentru norma 2, orice functie diferentiabila o data si tareconvexa satisface relatia:

f(x)f(y), x y x y22.


(iii) Functia g(x) = 12x22 este 1-tare convexa n raport cu norma 2.

(iv) Functia g(x) = 12x2p este (p 1)-tare convexa n raport cu norma

p, pentru p > 1.

(v) Functia g(x) = log nni=1

xi log xi este 1-tare convexa n raport cu

norma 1 pe domeniul {x Rn+ : x1 1}.Problema 16. Fie functia f definita de:

f(y, t) = maxxRn

y, x t2x2.

Determinati domeniul functiei f si aratati ca expresia acesteia se rescrien urmatoarea forma:

f(y, t) =

{0, daca y = 0, t = 0y22t, daca t > 0.

.

Problema 17. Sa se demonstreze ca pentru o matrice simetrica si pozitivdefinita Q Rnn si doi vectori u, v Rn are loc urmatoarea inegalitate:

|uTv| uTQuvTQ1v.

Problema 18. Fie matricele E,H si F de dimensiuni compatibile siF TF I. Sa demonstreze ca pentru orice > 0, avem

EFH +HTF TET EET + 1HTH.

Indiciu: Se foloseste inegalitatea (ET FH)T (ET FH) 0.

Problema 19. Fie matricea simetrica M de forma:

M =

[A BBT C

],

n care matricea C este inversabila. Se numeste complementul Schur allui M , matricea S = A BC1BT . Un rezultat des folosit n teoriamatriceala precizeaza ca urmatoarele relatii sunt echivalente:

(i) M 0 (M este pozitiv semidefinita);


(ii) A 0, (I AA)B = 0, C BTAB 0;(iii) C 0, (I CC)B = 0, A BTCB 0,unde cu A am notat pseudo-inversa matricei A. Mai mult, daca M 0obtinem un caz particular al echivalentelor:

(i) M 0 (M este pozitiv definita);(ii) A 0, C BTA1B 0;(iii) C 0, ABTC1B 0.Evidentiem, mai departe, un exemplu de aplicatie al acestui rezultat. Fiefunctia f : Rn Sn++ R definita de:

f(x, Y ) = xTY 1x.

Sa se demonstreze ca f este convexa folosind proprietatilecomplementului Schur.

Capitolul 3

Metode de ordinul I

3.1 Preliminarii

In acest capitol abordam probleme neliniare de optimizare neconstransa(unconstrained nonlinear programming - UNLP):

(UNLP ) : minxRn

f(x), (3.1)

unde functia obiectiv f este de doua ori diferentiabila. Conformconditiilor de optimalitate necesare pentru problema (3.1), orice punct deminim local pentru problema (UNLP), x, satisface urmatoarele relatii:

f(x) = 0, 2f(x)

64 Capitolul 3. Metode de ordinul I

I prezinta o complexitate a iteratiei foarte scazuta (n comparatie cucei de ordin II) si o convergenta accelerata n regiunile ndepartate depunctul de optim, atunci cand algoritmul intra n vecinatatea punctuluide optim, viteza acestora scade considerabil. De aceea, gasirea unuipunct de optim cu o acuratete mare este un proces dificil pentrumetodele de ordin I. In cazul problemelor de dimensiuni foarte mari,cand nu este necesara aflarea punctului de optim cu o acuratete ridicata,recomandarea principala pentru rezolvarea acestora sunt algoritmii deordin I datorita complexitatii reduse a iteratiilor acestora. In continuare,prezentam principalele metode de ordin I si exemple de functionare aleacestora.

3.2 Probleme rezolvate de laborator

3.2.1 Metoda Gradient

Metoda gradient se afla printre primele si cele mai simple metodedezvoltate n scopul determinarii unui punct critic aflat pe o anumitacurba (Cauchy, 1847). In principiu, metoda gradient reprezinta unalgoritm de ordin I care genereaza un sir de puncte (vectori) x1, x2, . . . ,pornind dintr-un punct initial ales. Structura esentiala a metodeigradient este enuntata n continuare:

Metoda Gradient

1. Se alege punctul initial x0, k := 0.

2. Se determina pasul k si se actualizeaza xk+1 = xk kf(xk).3. Daca criteriul de oprire nu este satisfacut, atunci se incrementeaza

k := k + 1 si se reia pasul 2,

undef(x) reprezinta gradientul functiei f n punctul x. Pentru alegereapasului k avem mai multe optiuni:

(i) Alegerea ideala a pasului k la fiecare iteratie presupune ca functiascalara () = f(xk f(xk)) sa descreasca cat mai mult posibil, i.e:

k = argmin>0

(),

numita si problema de line search.(ii) Deseori, n functie de f , minimizarea lui () poate fi foarte dificila.In acest caz, k poate fi gasit prin diversi algoritmi mai simpli de


cautare, ce includ conditii necesare asupra pasului pentru asigurarea uneidescresteri suficiente a functiei. Conditiile Wolfe reprezinta un exempluelocvent pentru aceasta strategie de line-search:

1. Se aleg doua constante c1 si c2 ce satisfac 0 < c1 < c2 < 1

2. Se determina k > 0 astfel ncat:

f(xk kf(xk)) f(xk) c1kf(xk)Tf(xk) (3.2)

f(xk)Tf(xk kf(xk)) c2f(xk)Tf(xk). (3.3)(iii) Un caz particular des utilizat n practica este metoda backtrackingce ajusteaza dimensiunea pasului k pentru ca prima relatie Wolfe (3.2)sa fie satisfacuta; metoda presupune alegerea unui parametru (0, 1]si actualizarea dimensiunii pasului, dupa cum urmeaza:

1. Se alege 0 > 0, (0, 1];2. Cat timp k nu satisface prima conditie Wolfe (3.2) iteram :

2.1. k+1 = k; k = k + 1.

(iv) Pentru functiile cu gradient continuu n sens Lipschitz cu o constantaL > 0, putem alege pasul k constant la fiecare iteratie. Daca aplicamiteratia metodei gradient, din conditia continuitatii functiilor cu gradientLipschitz avem:

f(xk+1) f(xk) k(1 L2k) f(xk)2 ,

rezultand ca trebuie sa selectam k (0, 2L), iar pentru o descrestereoptima a functiei, la fiecare iteratie alegem k =

1L.

In exemplul urmator vom implementa metoda gradient pentru prima sia treia dintre optiunile alegerii pasului k, unde criteriul de oprire va fiimpus de scaderea termenului f(xk) sub o precizie data.Exemplul 15. Fie functia f(x) = (x1 2)4 + (x1 2x2)2. Sa seimplementeze metoda gradient pentru rezolvarea problemei:

minxRn

f(x), (3.4)

n varianta cu pas ideal si cea cu pas ales prin metoda de backtracking.


Rezolvare. Pentru nceput, vom avea nevoie de doua functii:

f=feval_obj(x), g=gradient_obj(x)

care sa returneze valoarea functiei intr-un punct x, respectiv gradientulfunctiei n acel punct. Din moment ce vom cauta pasul ideal la fiecareiteratie a metodei gradient, va fi necesara o functie ce returneaza valoareafunctiei () = f(x+ d):

function f=phi_obj(alpha,x,d)

f=feval_obj(x+alpha*d);

end

Pentru gasirea pasului ideal la fiecare iteratie, vom utiliza functiafminsearch. Vom porni de la un punct initial x0, iar conditia deoprire a algoritmului va presupune ca norma gradientului sa fie sub oanumita toleranta impusa eps. Implementarea algoritmului este data deurmatoarea secventa de cod:

function xmin=gradient_method(x0,eps)

% Initializam vectori/matrice pentru

% memorarea gradientilor, a punctelor x

% generate de algoritm, etc

puncte_gradient=[]; puncte_iteratie=[];

valori_functie=[]; norme_gradienti=[];

%vom utiliza un vector g pentru a stoca gradientul curent

x=x0; g=gradient_obj(x);

while(norm(g)>eps)

g=gradient_obj(x); puncte_gradient=[puncte_gradient g];

puncte_iteratie=[puncte_iteratie x];

valori_functie=[valori_functie; feval_obj(x)];

norme_gradienti=[norme_gradienti; norm(g)];

alpha=fminsearch(@(alpha) phi_obj(alpha,x,-g).1);

x=x-alpha*g

end

xmin=x;

%Pentru afisarea grafica a rezultatelor,

%avem urmatoarele instructiuni

t=1:length(valori_functie);

figure(1)

hold on


plot(t,norme_gradienti(t),k,LineWidth,2);

hold off

figure(2)

hold on

plot(t,valori_functie(t),k,LineWidth,2);

hold off

%Pentru trasarea liniilor de contur si evolutia

%metodei gradient, avem urmatoarele

%instructiuni

[x1,x2]=meshgrid([1.2:0.01:2.8],[0.4:0.01:1.6]);

z=(x1-2).^4+(x1-2.*x2).^2;

figure(3)

hold on

contour(x1,x2,z,valori_functie);

plot3(puncte_iteratie(1,:),puncte_iteratie(2,:),...

valori_functie,r);

scatter3(puncte_iteratie(1,:),puncte_iteratie(2,:),...

valori_functie,filled);

hold off

1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.80.4

0.6

0.8

1

1.2

1.4

1.6

x1

x2

Figura 3.1: Convergenta metodei gradient cu pas ideal.

Apelarea functiei precedente se face n linia de comanda din Matlab, e.g.:

xmin=gradient_method([1.4;0.5],0.0001)

Pentru varianta metodei gradient cu pasul determinat de metoda debacktracking, se poate nlocui n cod functia fminsearch cu urmatoareasecventa de cod:


function alpha=backtrack_alpha(x,d)

alpha=1; t1=0.9; t2=0.2;

g=gradient_obj(x);

%Va trebui satisfacuta conditia Armijo:

while(feval_obj(x+alpha*d)>feval_obj(x)+t1*alpha*g*d)

alpha=alpha*t2;

end

care va fi apelata cu d = g = f(x).In Fig. 3.1 se observa caracteristica elocventa a metodei gradient carea fost precizata n sectiunile anterioare, si anume decelerarea ratei deconvergenta pe masura ce algoritmul se apropie de punctul de optim.In plus, Fig. 3.2 reda rezultatele grafice comparative ale convergenteimetodei gradient cand criteriul de oprire este de forma f(xk) f sauf(xk). Observam ca desi criteriul f(xk) f reflecta mult mai bineacuratetea punctului curent, acesta este rar folosit n practica, deoarecevaloarea optima nu se cunoaste a priori.

0 10 20 30 40 50 600

0.05

0.1

0.15

0.2

0.25

0.3

0.35

k

f(xk)f*

0 10 20 30 40 50 600

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

k

|| f(xk) ||

Figura 3.2: Comparatia convergentei variantelor metodei gradient (cucriteriul f(xk) f n prima figura si cu criteriul f(xk) n a doua),pentru pas ideal (linie continua) si pas obtinut prin backtracking (linie

punctata).

In exemplul urmator, vom implementa metoda gradient cu pas constantpentru o functie patratica cu gradient Lipschitz, pentru care L =max(Q).

Exemplul 16. Fie functia patratica:

f(x) =1

2xT

2.04 2.8 3.3;2.8 6 43.3 4 17.25

x+ [1 1 2]x


Figura 3.3: Progresul metodei gradient pornind dintr-un punct initial alesaleatoriu, spre punctul de optim.

Sa se implementeze metoda gradient pentru rezolvarea problemei:

minxR3

f(x).

Rezolvare. Implementarea metodei gradient din cerinta precedenta estedata de urmatoarea secventa de cod:

function [x,f,iter]=gradient(eps)

Q=[2.04 -2.8 3.3;...

-2.8 6 -4;...

3.3 -4 17.25];

q=[1;-1;-2];

L=max(eig(Q)); x=rand(3,1);

grad=Q*x+q;

f_new=0.5*x*Q*x+q*x; f_old=f_new+1; iter=0;

while ((f_old-f_new)>eps)

f_old=f_new; x=x-(1/L)*grad;

grad=Q*x+q; f_new=0.5*x*Q*x+q*x;

iter=iter+1;

end

f=f_new;

end

Sa se rezolve sistemul Qx = q si sa se compare solutia acestuia cu ceaa problemei anterioare. Sa se comenteze observatiile facute.


3.2.2 Metoda gradientilor con

Date post:	14-Oct-2015
Category:	Documents
Upload:	cristina-serban
View:	54 times
Download:	3 times

LabSem

Documents