Ett avståndsexempel

Följande exempel kommer också från en föreläsning. Problemtexten lyder: "Bestäm den punkt som ligger närmast origo och som samtidigt ligger på sfären, (x-1)^2+y^2+z^2=1, och i planet x+y+z=1."

Steg 1: karakterisera objektfunktion och bivillkor
På föreläsningen visade vi att det räcker att minimera x^2+y^2+z^2, kvadratroten behövs inte, med andra ord. Objektfunktionen är ickelinjär (en kvadratsumma, så vi kunde ha valt en specialrutin för denna typ av objektfunktion). Vi har två likhetsbivillkor, ett linjärt (planet) och ett ickelinjärt (sfären).

Steg 2: skriv objektfunktion och bivillkor på standardform
Objektfunktionen är redan på standardform, liksom det linjära bivillkoret. Vi skriver det ickelinjära villkoret som (x-1)^2+y^2+z^2-1=0 (optimeringsrutinen vill ha en nolla i ena ledet).

Steg 3: förpackning; vi inför de vektorer och matriser som behövs för att definiera objektfunktion och bivillkor
Vi inför först x som vektorn som lagrar variablerna, detta för att visa på ett potentiellt problem. x(1) svarar mot x (usch!), x(2) mot y och x(3) mot z. För att slippa göra fel byter vi beteckning och döper vektorn till w i stället, så w(1) svarar mot x, w(2) mot y och w(3) mot z

Likhetsbivillkoret definieras av radvektorn A_eq och skalären B_eq, där A_eq * w = Beq. Så, A_eq = [1, 1, 1]; och B_eq = 1; . Det ickelinjära bivillkoret definieras via en funktion, [in_eq, eq] = constr_fun(w), som hanterar både ickelinjära likhetsbivillkor och ickelinjära olikhetsbivillkor. Funktionen returnerar avvikelsen från bivillkoren, så i rutinen kommer vi att bilda

eq = (w(1)-1)^2 + w(2)^2 + w(3)^2 - 1;

Så, om avvikelsen, eq, är noll, så uppfyller w bivillkoret och om eq är skilt från noll så är bivillkoret inte uppfyllt.
Eftersom vi inte har något olikhetsbivillkor sätter vi in_eq till en tom variabel [] .

Steg 4: skriv Matlabkoden
Jag har valt att lägga de tre rutinerna i en fil, avst_ex.m .

function avst_ex
A    = [];               % inget linjärt olikhetsbivillkor
B    = [];

A_eq = [1, 1, 1];        % linjärt likhetsbivillkor
B_eq = 1;

LB   = [];               % inga enkla gränser
UB   = [];

w_guess = [1; 1; 0];     % start-gissning

[w_opt, obj_val] = fmincon(@obj_fun, w_guess, A, B, A_eq, B_eq, LB, UB, @constr_fun)

function obj_val = obj_fun(w)
obj_val = sum(w.^2);     % objektfunktion

function [in_eq, eq] = constr_fun(w)
in_eq = [];                                 % inget ickelinjärt olikhetsbivillkor
eq    = (w(1) - 1)^2 + w(2)^2 + w(3)^2 - 1; % ickelinjärt likhetsbivillkor

Körning (i redigerad form)

>> avst_ex
w_opt =
   1.8350e-01
   4.0825e-01
   4.0825e-01
obj_val =
   3.6701e-01

% facit från föreläsningen
>> s = 1 / sqrt(6);
>> [1-2*s, s, s]
ans =
   1.8350e-01   4.0825e-01   4.0825e-01  % stämmer

Antag nu att vi sätter w_guess = [1; -1; -1]; . Rutinen konvergerar fortfarande men hittar i stället ett lokalt maximum (det  största avståndet i själva verket; se föreläsningen). Vi har alltså inga garantier för att vi hittar ett lokalt minimum. Rutinen använder gradienten (bland annat) för att avgöra när den skall stanna, men gradienten ger oss inte information om vi har min eller max, för det krävs ju andraderivator (Hessian). Rutinen är en kvasi-Newton-rutin som approximerar (inversen av) Hessianen av Lagrangefunktionen.

>> avst_ex
w_opt =
   1.8165e+00
  -4.0812e-01
  -4.0838e-01
obj_val =
   3.6330e+00

Back