Lineair Programmeren in Mathematica

Simplexmethode

Luc Florack, 14 juni 2000

Stap 1: Specificeer de doelfunctie:

We definiëren de functie profit welke we willen maximaliseren:

winst[x_, y_, z_, u_, v_, w_] := 
      6360 + 6  x + 3  y + 8  z + 3  u + 9  v + 5  w

Stap 2: Specificeer het domein:

Het domein van deze functie is een polytoop, welke we definiëren middels een aantal ongelijkheden:

polytoop = {x ≥ 0, y ≥ 0, 
        z ≥ 0, u ≥ 0, v ≥ 0, 
        w ≥ 0, 480 - x - y ≥ 0, 
        400 - z - u ≥ 0, 230 - v - w ≥ 0, 
        x + z + v ≤ 420, y + u + w ≤ 250}

Stap 3: Optimaliseer:

Mathematica kent de functie ConstrainedMin. Aangezien deze minimaliseert moeten we -winst in plaats van winst nemen:

oplossing = 
      ConstrainedMin[-winst[x, y, z, u, v, w], 
        polytoop, {x, y, z, u, v, w}]

maximumwinst = -oplossing[[1]]

Alternatieven:

Mathematica kent trouwens ook de functie ConstrainedMax:

oplossing = 
      ConstrainedMax[winst[x, y, z, u, v, w], 
        polytoop, {x, y, z, u, v, w}]

maximumwinst = oplossing[[1]]

Een nog handiger methode is de volgende. Definiëer de lineaire vorm winstvorm als de coëfficiëntenvector van het lineaire deel van de doelfunctie; het minteken is wederom nodig omdat LinearProgramming minimaliseert:

winstvorm = -{6, 3, 8, 3, 9, 5}

Deze lineaire vorm is het eerste argument van  LinearProgramming. Als tweede argument voeren we een matrix in welke de coëfficiënten bevat van de ongelijkheden die de polytoop definiëren. Ook hier moeten we op de tekens letten; Mathematica gaat standaard uit van "groter dan" ongelijkheden. Als derde argument geven we een rijtjesvector met de constanten uit de ongelijkheden, wederom voorzien van het juiste teken:

polytoopmatrix=(({{-1, -1, 0, 0, 0, 0}, {0, 0, -1, -1, 0, 0}, {0, 0, 0, 0, -1, -1}, {-1, 0, -1, 0, -1, 0}, {0, -1, 0, -1, 0, -1}}))

polytoopvector = {-480, -400, -230, -420, -250}

oplossing = 
      LinearProgramming[winstvorm, polytoopmatrix, polytoopvector]

Vergeet bij het bepalen van de maximale winst niet de offset uit de doelfunctie, en let weer op het teken:

winstoffset = 6360

maximumwinst = -winstvorm . oplossing  + winstoffset

Tips:

Met behulp van de functie UnitStep kun je een polytoop van willekeurige dimensie definiëren als een binaire functie, en zodoende gemakkelijk testen of een bepaald punt aanvaardbaar is of niet. In 2 dimensies kun je polytopen op deze wijze tevens visualiseren, al moet je wel veel punten plotten om een enigszins fatsoenlijk plaatje te krijgen:

polyeder[x_, y_] := 
      UnitStep[x] UnitStep[y] UnitStep[-x + y - 1] UnitStep[x + y - 3] 
        UnitStep[-2 x - y + 4]

{polyeder[1, 4] == 1, polyeder[0, 4] == 1}

DensityPlot[polyeder[x, y], {x, -0.5, 2}, {y, 1.5, 5}, 
        PlotRange -> All, Frame -> True, PlotPoints -> 250, 
        Mesh -> False, DefaultFont -> {"Helvetica-Bold", 16. }, 
        PlotLabel -> ""];