Paketexemplet

I detta exempel attackerar vi paketproblemet från föreläsningen. Det är dock två väsentliga skillnader. Under föreläsningen bestämde vi, med handräkning, den största volymen. I detta exempel använder vi fmincon för att bestämma ett lokalt maximum (som råkar vara globalt maximum).

Vi vill alltså skicka ett paket, med maximal volym, till utlandet (exemplet gjordes 2005 och villkoren kan ha ändrats). I postens ``måttguide för paket i kartong, utland'' fanns följande olikhetsbivillkor angivna, om l är längd, b bredd och h höjd:

l <= 150 cm
l+2(b+h) <= 300 cm

Vår objektfunktion är f(l,b,h) = lbh, som skall maximeras. För att få ett kompakt område lade vi till tre villkor:
0 <= l, b, h

För att lösa problemet med Matlab utför vi fyra steg:

Steg 1: karakterisera objektfunktion och bivillkor
Objektfunktionen är ickelinjär (hade den varit linjär hade vi i stället valt rutinen linprog som är speciellt anpassad för linjär objektfunktion och linjära bivillkor). Alla bivillkoren är linjära olikhetsbivillkor och alla utom l+2(b+h) <= 300 cm är enkla gränser.

Steg 2: skriv objektfunktion och bivillkor på standardform
Optimeringsrutiner brukar minimera objektfunktionen, men vi vill ju maximera densamma. Detta åstadkommer vi genom att minimera -lbh. Nu till bivillkoren. De enkla gränserna kan skrivas:

0 <= l
0 <= b
0 <= h

l <= 150
b <= Inf   b eller h kan givetvis inte bli oändligt stora, men vi har inga
h <= Inf   uppenbara övre gränser för bredd eller höjd (*)

l+2(b+h) <= 300 skriver vi om som  l+2b+2h <= 300

(*) Jag skrev detta mest för att visa hur man kan använda Inf (och -Inf för en undre gräns). I detta fall är det enkelt att härleda övre gränser med hjälp av l+2(b+h) <= 300. b och h kan maximalt vara 150 (om övriga variabler är noll). I labben går det inte bra med -Inf, Inf. Du måste ge rutinen mer hjälp (information).

Steg 3: förpackning; vi inför de vektorer och matriser som behövs för att definiera objektfunktion och bivillkor
Optimeringsrutinen arbetar med en vektor av obekanta (och inte med skalära värden l, b, h). Vi förpackar därför l, b, h i en vektor, kalla den x (men den kan ju givetvis heta något annat). Vi bestämmer oss också för en ordning, x(1) svarar mot l, x(2) mot b och x(3) mot h. Det går bra med en annan ordning, men man måste givetvis komma ihåg vilket element i x som svarar mot en given variabel. Nu till bivillkoren.

För de undre, enkla gränserna använder vi kolonnvektorn LB (lower bound) som definieras som LB = [0; 0; 0]; och för de övre enkla gränserna har vi kolonnvektorn UB (upper bound) som definieras som LB = [150; Inf; Inf]; Man kan givetvis skapa vektorerna på andra sätt. LB = zeros(3, 1); eller LB = [0, 0, 0]'; är två andra alternativ. För det sista bivillkoret bryter vi ut koefficienterna och placerar dessa i en radvektor, A (för att använda Matlabs variabelnamn). Högerledet, 300, placerar vi i variabeln B. Bivillkoret kan då skrivas A * x <= B, där * står för vanlig matrismultiplikation. Så, A = [1, 2, 2]; och B = 300; .

Steg 4: skriv (och kör) Matlabkoden
Jag har valt att skriva en separat funktion för objektfunktionen. Detta behövs inte i enkla fall, men i labben är det praktiskt.


Låt oss lägga huvudprogrammet i paket.m .

LB   = [  0;   0;   0];  % enkla gränser
UB   = [150; Inf; Inf];

A    = [1, 2, 2];        % linjärt olikhetsbivillkor
B    = 300;

A_eq = [];               % inget linjärt likhetsbivillkor
B_eq = [];

x_guess = [10; 10; 10];  % start-gissning

% Anropa fmincon, x_opt är lösningen (om det går bra)
% och obj_val är objektfunktionens värde
[x_opt, obj_val] = fmincon(@obj_fun, x_guess, A, B, A_eq, B_eq, LB, UB)

Här följer objektfunktionen i en egen fil, obj_fun.m

function obj_val = obj_fun(x)
obj_val = -x(1) * x(2) * x(3);   % OBS: minus. obj_val = -prod(x); går också bra

Körningen ser ut väsentligen ut så här (jag har tagit bort några ointressanta rader):

>> paket
Optimization terminated: magnitude of directional derivative in search
 direction less than 2*options.TolFun and maximum constraint violation
  is less than options.TolCon.

Active inequalities (to within options.TolCon = 1e-06):
  lower      upper     ineqlin   ineqnonlin
                          1
x_opt =
   1.0000e+02
   5.0000e+01
   5.0000e+01
obj_val =
  -2.5000e+05

Låt oss jämföra detta med föreläsningen. Det står:
"Vi jämför nu  de inramade kvantiteterna och finner att maximum är 250000 vilket antas för (l,b,h)=(100,50,50)." I utskriften ovan framgår att fmincon har hittat samma värden (vi får givetvis ta bort minustecknet från obj_val). Avsnittet med "Active inequalities" säger att minimipunkten ligger på randen av ett av olikhetsbivillkoren (ineqlin = 1), dvs. punkten ligger på planet som definieras av l+2(b+h) = 300, detta konstaterade vi också på föreläsningen. Olikhetsbivillkoret är alltså uppfyllt med likhet (det finns ingen tabell med likhetsbivillkor, eftersom dessa alltid skall vara uppfyllda med likhet). Raden om "terminated" talar om att rutinen konvergerade och hittade en punkt där konvergenskriteriet är uppfyllt.

Vi kan paketera huvudprogram och objektfunktion i en enda fil, förutsatt att vi gör om huvudprogrammet till en funktion och att denna kommer först i filen. Så, i filen paket.m lägger vi

function paket           % OBS: nu en funktion
LB   = [  0;   0;   0];  % enkla gränser
UB   = [150; Inf; Inf];

A    = [1, 2, 2];        % linjärt olikhetsbivillkor
B    = 300;

A_eq = [];               % inget linjärt likhetsbivillkor
B_eq = [];

x_guess = [10; 10; 10];  % start-gissning

[x_opt, obj_val] = fmincon(@obj_fun, x_guess, A, B, A_eq, B_eq, LB, UB)

% Här följer koden för obj_fun i samma fil
function obj_val = obj_fun(x)
obj_val = -x(1) * x(2) * x(3);   % OBS: minus. obj_val = -prod(x); går också bra

Back