Laboratorium Wydział Elektroniki	Technika optymalizacji
	Data wykonania ćwiczenia: 5.12.2012 Prowadząca: Dr inż. Ewa Szlachcic
Ćwiczenie IV Zadanie programowania nieliniowego z ograniczeniami.

1. Wstęp teoretyczny – opis metody kar i barier.

Metoda gradientowa służy do rozwiązywania zadań programowania nieliniowego. Za jej pomocą rozwiązuje się zadania, w których chcemy odszukać najmniejszą wartość funkcji celu. Kierunek poszukiwań (poprawy) to kierunek, w którym maleje wartość funkcji celu. W punkcie x^(k) generuje się kierunek poprawy d^(k), spełniający warunek: .

Najczęściej jako kierunek poszukiwań rozwiązania minimalnego wybiera się , tzn. kierunek największego spadku funkcji w punktcie x^(k). Dla każdego wznaczonego punktu sprawdza się czy spełnia on zadane ograniczenia. Kryterium zatrzymania algorytmu to: dla ograniczeń nieaktywnych, lub dla ograniczeń aktywnych.

1. Rozwiązania zadań w Matlabie.

f(x) = x₁ ² + x₁ * x₂ + 0.5 * x₂ ² − x₁ − x₂, punkt startowy x₀ = [3,  3]

Punkt optymalny: $\hat{x} = \left\lbrack 0,\ 1 \right\rbrack,$ $f(\hat{x}) = - 0.5$

1. Ograniczenia 1 : x₁−x₂≥1

x₁+x₂≤4

Kod z Matlaba

[x y]=meshgrid(-4:.1:4);

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

con =-x+y+1;

con1 = x+y-4;

contour(x,y,z)

hold on

contour(x,y,con,[1e-8])

hold on

contour(x,y,con1,[1e-8])

Kontur funkcji celu i obszar:

Rozwiązanie za pomocą funkcji fmincon dla tolerancji tolcon=10^-1, tolfun=10^-1, punkt startowy to [2 -2], znajdującu się wewnątrz obszaru ograniczeń.

Kod z Matlaba:

tolfun=10^-1;

tolcon=10^-1;

options=optimset('TolFun', tolfun, 'TolCon',tolcon,'PlotFcns','optimplotfval', 'display' , 'iter');

fun = @(x)x(1).^2+x(1).*x(2)+0.5.*x(2).^2-x(1)-x(2);

x0=[2 -2];

A=[-1 1; 1 1];

b=[-1 ;4];

[x,fval,exitflag,output]=fmincon(fun,x0,A,b,[],[],[],[],[],options)

Wyniki z Matlaba

Max Line search Directional First-order

Iter F-count f(x) constraint steplength derivative optimality Procedure

0 3 2 -3

1 6 0.5 -1 1 -1.41 1

2 9 -2.48353e-009 0 1 -1 0.5

3 12 -0.1 0 1 -0.707 5.12e-009

Optimization terminated: first-order optimality measure less

than options.TolFun and maximum constraint violation is less

than options.TolCon.

Active inequalities (to within options.TolCon = 0.1):

lower upper ineqlin ineqnonlin

1

x =

0.8000 -0.2000

fval =

-0.1000

exitflag =

1

output =

iterations: 3

funcCount: 12

lssteplength: 1

stepsize: 0.2828

algorithm: 'medium-scale: SQP, Quasi-Newton, line-search'

firstorderopt: 5.1161e-009

constrviolation: 0

message: [1x144 char]

Sprawdzenie czy wyliczony punkt mieści się w ograniczonym obszarze:

0,8-(-0,2)≥1 PRAWDA

0,8+(-0,2)≤4 PRAWDA

Wyliczony punkt mieści się w obszarze ograniczeń.

Metodą prób i błędów, doszliśmy do wniosku, iż najlepszym punktem startowym będzie punkt [2 -2] . Obecność kolumny max constraint w wyniku z MATLABa pokazuje, że zadany przykład obliczeniowy zawiera ograniczenia i algorytm w każdym kroku iteracji sprawdza czy nowo uzyskany punkt należy do dozwolonego obszaru rozwiązań. Jeżeli tak jest, to wartość max constraint jest <=0, w przeciwnym wypaku znak zmienia się na dodatni. Przy zbliżaniu się do ograniczenia, każda kolejna iteracja zmniejsza funkcję o coraz mniejszą wartośc. Minimalne oddalenie osiągamy przy 3 iteracji .Punkt znajduje się w obszarze rozwiązań, jest równy zero, a wartość funkcji wynosi -0,1. W przypadku tego przykładu, dla każdej wartości tolfun oraz tolcon mniejszej lub równej 10^-1 otrzymywaliśmy ten sam wynik. Jeśli wartośc tą ustawiliśmy na wiekszą od 1, wynik był znacznie oddalony od optymalnego rozwiązania.

1. Ograniczenie 2 : x₂≤x₁²

rozpatrywane są każde punkty znajdujące się poza wnętrzem paraboli

Jako punkt startowy wybraliśmy [0.5 0.2], tolerancja tolcon=10^-1

Kod z Matlaba:

[x y]=meshgrid(-4:.1:4);

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

con =-x+y+1;

con1 = x+y-4;

con2=y-x.^2;

contour(x,y,z)

hold on

contour(x,y,con2,[1e-8])

tolcon=10^-1;

options=optimset( 'TolCon',tolcon,'PlotFcns',@optimplotconstrviolation, 'display' , 'iter');

fun = @(x)x(1).^2+x(1).*x(2)+0.5.*x(2).^2-x(1)-x(2);

x0=[0.5 0.2];

[x,fval,exitlog,out] = fmincon(fun,x0,[],[],[],[],[],[],@(x) mycon(x),options)

Funkcja mycon.m

function [c,ceq] = mycon(x)

c=x(2)-x(1)^2;

ceq=[];

Wyniki z Matlaba:

Max Line search Directional First-order

Iter F-count f(x) constraint steplength derivative optimality Procedure

0 3 4.5 -12

1 6 1.91417 -3.378 1 -2.24 1.84

2 9 -0.269996 -0.2526 1 -1.66 0.219

3 12 -0.338515 -0.01789 1 -0.514 0.0684

4 15 -0.343592 -0.0006294 1 -0.225 0.00146

5 18 -0.34375 -6.188e-007 1 -0.198 2.31e-006

Optimization terminated: first-order optimality measure less

than options.TolFun and maximum constraint violation is less

than options.TolCon.

Active inequalities (to within options.TolCon = 0.1):

lower upper ineqlin ineqnonlin

1

x =

0.5000 0.2500

fval =

-0.3437

exitlog =

1

out =

iterations: 6

funcCount: 18

lssteplength: 1

stepsize: 6.3861e-007

algorithm: 'medium-scale: SQP, Quasi-Newton, line-search'

firstorderopt: 6.6916e-007

constrviolation: -6.1876e-007

Sprawdzenie czy obliczony punkt optymalny znajduje się w rozpatrywanym obszarze:

2500≤0, 5000² PRAWDA

Obliczony punkt znajduje się na linii oganiczającej, tak więc należy do analizowanego obszaru i spełnia ograniczenia aktywne z dokładnością 10^-4

Jako punkt startowy wybraliśmy [-10 -10], dokłaność tolcon=10^-1

Wynik Matlaba:

Max Line search Directional First-order

Iter F-count f(x) constraint steplength derivative optimality Procedure

0 3 270 -110

1 6 19.6305 -20.4 1 -27.4 27

2 9 5.81759 -4.947 1 -3.9 1.52

3 12 1.65886 -1.555 1 -1.99 0.949

4 15 0.771908 -0.3169 1 -1.16 0.423

5 18 0.547118 -0.04345 1 -0.703 0.098 Hessian modified

6 21 0.509355 -0.0047 1 -0.29 0.0515 Hessian modified

7 24 0.504381 -0.0002433 1 -0.147 0.303 Hessian modified

8 27 0.504288 -0.003752 1 -0.0337 0.262 Hessian modified

9 31 0.503278 -0.003268 0.5 -0.0296 0.0401

10 34 0.50163 -0.001653 1 -0.0364 0.0271

11 37 0.500053 -5.193e-005 1 -0.109 0.0236 Hessian modified

12 40 0.500001 -1.26e-007 1 -0.0682 0.000315 Hessian modified

Optimization terminated: directional derivative predicts change in

objective value less than options.TolFun and maximum constraint violation

is less than options.TolCon.

Active inequalities (to within options.TolCon = 0.1):

lower upper ineqlin ineqnonlin

1

x =

-1.0103 1.0208

fval =

0.5000

exitlog =

5

out =

iterations: 13

funcCount: 40

lssteplength: 1

stepsize: 0.0014

algorithm: 'medium-scale: SQP, Quasi-Newton, line-search'

firstorderopt: 1.2789e-004

constrviolation: -1.2601e-007

message: [1x173 char]

Sprawdzenie czy obliczony punkt optymalny znajduje się w rozpatrywanym obszarze:

1, 0208 ≤ 1, 0207 FAŁSZ

Obliczony punkt lezy w obszarze ograniczeń. Prawdopodobnie wynika to z powodu wyboru zbyt odległego punktu startowego

W przypadku ograniczenia parabolą, obliczenie optymalnego punktu jest trudniejsze I mniej dokładne. Jako punkt startowy wybraliśmy oszacowane na podstawie warstwic minimum. Dla dokładności tolcon=10^-1 , program wykonał siedem iteracji, a otrzymany wynik był satysfakcjonujący. Każda próba zmniejszenia parametru tolcon dawała dokładnie ten sam wynik. Oznacza to, że dla funkcji wypyklych tolerancja spełnienia ograniczeń nie ma większego wplywu na dokładność obliczeń. Znaczne oddalenie punktu startowego od szacowanego minimum, powodowało wyliczenie innego punktu minimalnego, który mimo zadanej małej tolerancji, nie spełniał nierówności ograniczenia.

3.1. Przykład 17 z ograniczeniami

3.2. Ograniczenie 1

x₁ + x₂ ≥ 6

Kod z Matlaba:

[x y]=meshgrid(0:.1:6,0:.1:6);

z=4*x.^2-2.1*x.^4+1/3*x.^6+x.*y-4*y.^2+4*y.^4;

con=-x-y+6

figure(2)

contour(x,y,z)

hold on

contour(x,y,con,[1e-8])

x0=[3.2 3.2];

x0=[3.2 3.2];

tolcon=10^-1;

options=optimset('TolCon',tolcon,'PlotFcns',@optimplotfval, 'display' , 'iter');

fun = @(x)4.*x(1).^2-2.1.*x(1).^4+(1/3).*x(1).^6+x(1).*x(2)-4.*x(2).^2+4.*x(2).^4;

[x,fval,exitlog,out] = fmincon(fun,x0,[],[],[],[],[],[],@(x) mycon1(x),options)

Funkcja mycon1.m

function [c,ceq] = mycon1(x)

c=-x(1)-x(2)+6;

ceq=[];

Wyniki z Matlaba dla punktu startowego [3.2 3.2] oraz tolerancj tolcon=10^-1

Max Line search Directional First-order

Iter F-count f(x) constraint steplength derivative optimality Procedure

0 3 567.383 -0.4

1 15 564.887 -0.3992 0.00195 -58 185

2 18 398.421 -3.13e-009 1 -659 114

3 23 398.412 -2.348e-009 0.25 -3.13 1.05

4 26 398.412 2.631e-012 1 -0.243 0.000668

Optimization terminated: magnitude of search direction less than 2*options.TolX

and maximum constraint violation is less than options.TolCon.

Active inequalities (to within options.TolCon = 0.1):

lower upper ineqlin ineqnonlin

1

x =

3.1184 2.8816

fval =

398.4121

exitlog =

4

out =

iterations: 5

funcCount: 26

lssteplength: 1

stepsize: 6.2823e-007

algorithm: 'medium-scale: SQP, Quasi-Newton, line-search'

firstorderopt: 2.4372e-004

constrviolation: 2.6308e-012

message: [1x142 char]

Sprawdzenie czy obliczony punkt spełnia nierówność ograniczającą:

3.1184 +2.8816≤6

3.1184 +2.8816=6 PRAWDA

Obliczony punkt znajduje się na linii nierówności ograniczającej.

W przypadku ograniczeń liniowych, tolerancja spełnienia ograniczeń nie ma wplywu na wynik obliczeń. Zarówno dla tolcon=10^-1 jak i tolcon=10^-6 program wykonał dokładnie tyle samo iteracji, a wartość funkcji w poszczególnych krokach była taka sama w obu przypadkach.

3.2. Ograniczenie 2

Ograniczenie: (x₁−5)²+(x₂−5)²≤1

Rozpatrywanym obszarem jest wnętrze okręgu o promieniu 1

Program z Matlaba

[x y]=meshgrid(4:.1:6,4:.1:6);

z=4*x.^2-2.1*x.^4+1/3*x.^6+x.*y-4*y.^2+4*y.^4;

con=(x-5).^2+(y-5).^2-1;

figure(2)

contour(x,y,z)

hold on

contour(x,y,con,[1e-8])

x0=[5 5];

tolcon=10^-1;

options=optimset('TolCon',tolcon,'PlotFcns',@optimplotconstrviolation, 'display' , 'iter');

fun = @(x)4.*x(1).^2-2.1.*x(1).^4+(1/3).*x(1).^6+x(1).*x(2)-4.*x(2).^2+4.*x(2).^4;

[x,fval,exitlog,out] = fmincon(fun,x0,[],[],[],[],[],[],@(x) mycon2(x),options)

Funkcja mycon2.m

function [c,ceq] = mycon2(x)

c=(x(1)-5)^2+(x(2)-5)^2-1;

ceq=[];

Wyniki z Matlaba dla punktu startowego [5 5] oraz tolerancji tolcon=10^-1

Optimization terminated: magnitude of search direction less than 2*options.TolX

and maximum constraint violation is less than options.TolCon.

Active inequalities (to within options.TolCon = 0.1):

lower upper ineqlin ineqnonlin

1

x =

4.1773 4.4315

fval =

2.6841e+003

exitlog =

4

out =

iterations: 20

funcCount: 94

lssteplength: 1

stepsize: 1.3987e-006

algorithm: 'medium-scale: SQP, Quasi-Newton, line-search'

firstorderopt: 0.0043

constrviolation: -1.3498e-006

message: [1x142 char]

Sprawdzenie czy otrzymany punkt spełnia ograniczenia za pomoca program Matlab

(4, 1773 − 5)²+(4, 4315 − 5)²≤1

>> (x(1)-5)^2+(x(2)-5)^2

ans =

1.0000 PRAWDA!

Otrzymany punkt mieści się w zadanym ograniczeniu.

Wyniki z Matlaba dla punktu startowego [5 5] oraz tolerancji tolcon=10^-8, były identyczne jak dla tolerancji tolcon = 10^-1. Z tego względu postanowiliśmy zrezygnowac z zamieszczania ich w tym sprawozdaniu.

4. Wnioski

W ćwiczeniu tym, przede wszystkim analizowaliśmy w programie MATLAB wpływ ograniczeń na poprawność rozwiązania. Zastosowany algorytm dobiera kolejne punkty na podstawie gradientu badanej funkcji w punkcie. Poruszając się w kierunku minimaizuje wartość funkcji.

W każdym kroku sprawdza czy nowo uzyskany punkt należy do dozwolonego obszaru rozwiązań. Jeżeli tak jest, to wartość max constraint jest ujemna, w przeciwnym wypaku znak zmienia się na dodatni. Wówczas algorytm dobiera kolejne kroki w taki sposób aby minimalizować wartość max constraint(C).

Warunkiem stopu jest lub . Spełnienie pierwszego z tych warunków oznacza, że algorytm wyznaczył minimum lokalne i ograniczenia były nieaktywne, natomiast spełnienie drugiego świadczy o aktywności ograniczeń i wyznaczeniu punkt, dla którego funkcja celu przyjmuje wartość minimalną na danym ograniczonym obszarze.