Simplexmethode
De simplexmethode (of het simplexalgoritme) is een methode in de wiskundige optimalisatie (zie ook operationeel onderzoek). De techniek werd in 1947 door George Dantzig ontwikkeld. De simplexmethode lost een lineaire-optimaliseringsprobleem in een eindig aantal stappen op, of stelt de onoplosbaarheid van het probleem vast. In theoretische gevallen kunnen cycli optreden, die het vinden van de optimale oplossing verhinderen. De naam komt van het feit dat de vergelijkingen van het probleem een simplex beschrijven, waarvan de rand bij het vinden van de oplossing beschreven wordt.
Geschiedenis
[bewerken | brontekst bewerken]De Russische wiskundige Leonid Kantorovich behandelde in 1939 als eerste lineaire optimalisatie in zijn boek "Wiskundige Methoden in de Organisatie en Planning van de Productie". Kort daarna presenteerde de Amerikaan F.L. Hitchcock een werk over een transportprobleem. Toen erkende men echter nog niet de betekenis van dit werk.
De doorbraak voor de simplexmethode kwam er door George Dantzig in 1947, toen ook hij een werk over lineaire optimalisatie publiceerde. In eerste instantie toonden het Amerikaanse Leger, in het bijzonder de US Air Force, interesse in dit werk - men wilde de militaire acties optimaliseren. De daaropvolgende jaren ontwikkelden John von Neumann en Oskar Morgenstern het algoritme verder.
Met de opkomst van computers halverwege de jaren 50 kon men grotere problemen oplossen. Men ontwikkelde speciale varianten van de simplexmethode die op de eigenschappen van de computers uit die tijd afgesteld waren. Ze gingen bijvoorbeeld heel zuinig om met het gebruik van het hoofdgeheugen.
Gelijktijdig - na 1950 - ontdekte de industrie dat men met de simplexmethode ook optimalisatieproblemen in hun domein kon oplossen. In het bijzonder gebruiken olieraffinaderijen de methode. In 1960 optimaliseerde W. Knödel de toelevering van de suikerfabrieken in Oostenrijk. Daarmee verlaagde hij de transportkosten met ongeveer 10 procent.
Definitie van het lineair optimalisatieprobleem
[bewerken | brontekst bewerken]- Maximaliseer de doelfunctie
- onder de nevenvoorwaarden en .
- Daarbij zijn
- A: m×n-matrix
- a: n-dimensionale vector
- b: m-dimensionale vector
- gegeven parameter, en de doelfunctie
- is als scalair product een lineaire functie in .
- Daarbij zijn
De doelfunctie wordt ook wel objectief genoemd, de nevenvoorwaarden ook wel restricties.
is hier de notatie van een matrixvermenigvuldiging. De 'ongelijkheden' en moeten component per component opgevat worden, bijvoorbeeld , met , betekent hetzelfde als .
Die punten (vectoren) die aan alle nevenvoorwaarden tegelijkertijd voldoen (dit zijn geldige oplossingen van het probleem) vormen een convexe polyeder in de n-dimensionale ruimte, een zogeheten simplex. Wanneer er geen geldige oplossingen zijn, spreken de nevenvoorwaarden elkaar tegen en is het probleem niet oplosbaar.
Wanneer de verzameling geldige oplossingen niet leeg is, zijn er twee mogelijkheden:
- De doelfunctie G heeft in de toegelaten punten een bovengrens:
Hier kan men zien dat het maximum van de doelfunctie door een punt op de rand van de geldige verzameling ingenomen wordt, en wel door een hoekpunt. Het simplexalgoritme verplaatst zich, vertrekkend van een geldige beginoplossing, langs de hoekpunten van de polyeder naar zijn optimale oplossing. - De doelfunctie G is niet naar boven begrensd:
In dat geval stelt men vast dat er een oneindig lange zijde van de polyeder met oneindige toenemende doelfunctie bestaat. Men zegt dat de optimale waarde van de doelfunctie is, en men kan geen optimale oplossing daarvoor geven.
Wanneer geen geldige beginoplossing gekend is, kan men een lichtjes aangepast lineair optimaliseringsprobleem gebruiken, om een beginoplossing te vinden, en voor de vaststelling van de oplosbaarheid van het probleem. (De zogenoemde Fase I met kunstmatige variabelen).
Bij veel probleemstellingen komt het voor dat men de doelfunctie G niet wil maximaliseren, maar wil minimaliseren. Bijvoorbeeld moeten bij een transportprobleem de transportkosten zo klein mogelijk zijn, of moet de wegen zo kort mogelijk zijn. In dit geval voert men eenvoudigweg een nieuwe doelfunctie in, en men maximaliseert deze.
Voorbeeld
[bewerken | brontekst bewerken]Aan de hand van een eenvoudig voorbeeld wordt de oplossingsmethode stap voor stap getoond. Reële problemen kunnen meer dan duizend randvoorwaarden hebben, waar men niet onmiddellijk eventuele tegenstellingen in kan herkennen, of dus niet direct het bestaan van een geldige beginoplossing kan vinden.
Probleemstelling in woorden
[bewerken | brontekst bewerken]Een firma biedt 2 verschillende producten aan. Er staan 3 machines A, B en C ter beschikking. Machine A heeft een maximale maandelijkse looptijd (capaciteit) van 170 uur, machine B van 150 uur en machine C van 180 uur. Een eenheid van product 1 geeft een winstbijdrage van 300 euro, een eenheid van product 2 daarentegen 500 euro. De vaste kosten bedragen 36.000 euro per maand.
Als men 1 eenheid van product 1 vervaardigt, dan gebruikt men daarvoor 1 uur machine A en daarna 1 uur machine B. 1 eenheid van product 2 vereist achtereenvolgens 2 uur machine A, 1 uur machine B en 3 uur machine C.
Wiskundige formulering met ongelijkheden
[bewerken | brontekst bewerken]Stel dat het bedrijf per maand x1 eenheden van product 1 en x2 eenheden van product 2 vervaardigt, dan bedraagt de opbrengst (winst of verlies)
G = 300x1 500x2 - 36.000
Deze opbrengst wil het bedrijf maximaliseren.
Aangezien de capaciteit van de machines begrensd is, kan men geen willekeurig grote hoeveelheden produceren. Uit de bezettingstijden leidt men de volgende nevenvoorwaarden voor x1 en x2 af:
1•x1 2•x2 ≤ 170 machine A 1•x1 1•x2 ≤ 150 machine B 0•x1 3•x2 ≤ 180 machine C
Negatieve productiehoeveelheden zijn niet mogelijk, dus gelden nog de nevenvoorwaarden
x1 ≥ 0 x2 ≥ 0
Omzetten in gelijkheden
[bewerken | brontekst bewerken]Omdat het met stelsels gelijkheden gemakkelijker rekenen is, worden de ongelijkheden in gelijkheden omgezet. Daartoe voert men de zogenoemde slackvariabelen yA, yB en yC in, die de niet benutte tijd van de respectievelijke machines voorstellen. Men vormt bijvoorbeeld de nevenvoorwaarde voor machine A in volgende gelijkheid om
yA x1 2•x2 = 170
Vanzelfsprekend mogen de slackvariabelen niet negatief zijn (niet-negativiteitsvoorwaarde).
Het probleem kan nu op volgende manier geformuleerd worden:
Maximaliseer de doelfunctie G onder de volgende nevenvoorwaarden:
G - 300x1 - 500x2 = - 36.000 yA x1 2x2 = 170 yB x1 x2 = 150 yC 3x2 = 180
yA, yB, yC, x1, x2 ≥ 0
Het Simplex-Tableau
[bewerken | brontekst bewerken]Deze vergelijkingen plaatst men in een zogeheten simplex-tableau, dat een verkorte voorstelling van het bijhorende stelsel van vergelijkingen is. Alle volgende bewerkingen zijn wiskundige berekeningen met het stelsel vergelijkingen.
-------------- ------------- | x1 x2 | rechterlid ---- -------------- ------------- G | -300 -500 | -36000 ---- -------------- ------------- YA | 1 2 | 170 = b1 YB | 1 1 | 150 = b2 YC | 0 3 | 180 = b3
De variabelen in de titelrij (x1, x2) heten niet-basisvariabelen, de variabelen in de eerste kolom basisvariabelen (YA, YB, YC). De getallen in "rij G" - de vergelijking van de doelfunctie - zijn de coëfficiënten van de doelfunctie. De variabelen b1, b2 en b3 bevatten de waarde van de rechterleden.
Bepaling van een geldige beginoplossing
[bewerken | brontekst bewerken]Vervolgens bepaalt men een geldige oplossing. Dat zijn waarden voor x1 en x2 die aan alle nevenvoorwaarden voldoen.
Men kan in dit geval onmiddellijk een triviale oplossing geven, namelijk
x1 = 0 x2 = 0
Deze oplossing wordt in het bovenstaande tableau opgenomen. x1 en x2 zijn de niet-basisvariabelen en hebben steeds de waarde "0". De basisvariabelen hebben de waarde van de onbenutte machinelooptijden. Aangezien met de beginoplossing niets geproduceerd wordt: YA = 170, YB = 150, YC = 180. Deze beginoplossing is natuurlijk onbevredigend: er wordt niet geproduceerd, de opbrengst G bedraagt -36.000 euro, wat verlies betekent. Daarom moet men proberen om betere geldige oplossingen te vinden.
Verbetering van de oplossing door middel van een simplex-iteratie
[bewerken | brontekst bewerken]In een iteratie van het simplexalgoritme wordt een basisvariabele met een niet-basisvariabele omgewisseld. Dit is een verwisselingsstap.
In aanmerking komen niet-basisvariabelen die een negatieve coëfficiënt in de doelfunctie hebben. Tussen deze variabelen zoekt men die variabele die bij omwisseling de grootste toename van de doelfunctie levert. Stel k het nummer van de kolom van de te verwisselen niet-basisvariabele, en r het rijnummer van de te verwisselen basisvariabele. De verwisselingsstap komt overeen met juist één stap bij het oplossen van een stelsel van gelijkheden, waarbij men rij r oplost naar de variabelen xk en dan xs in de rest van de gelijkheden inzet. Stel aij een matrix-element uit de simplex-tableau, dan heet ark het Pivotelement van de tableau. Kolom k is de pivotkolom, rij r de pivotrij.
In de verwisselingsstap berekent men het nieuwe tableau op volgende wijze:
Pivotelement:
- (Formule 1)
Pivotrij voor j k:
- (2)
- (3)
Pivotkolom voor i r:
- (4)
overige elementen:
- (5)
- (6)
Keuze van de pivotelementen
[bewerken | brontekst bewerken]Formule (6) is bepalend voor de keuze van de pivotrij r en de pivotkolom k. Men zoekt die waarden voor r en k waarvoor de coëfficiënt van de doelfunctie met index k negatief is, en waarvan de waarde van
zo groot mogelijk wordt. Daarbij dient men op te merken dat na de verwisseling de ontstane oplossing nog steeds geldig moet zijn, dit betekent dat alle nevenvoorwaarden vervuld moeten blijven. Daartoe kiest men r en k zodanig dat het quotiënt
- met
binnen een kolom minimaal wordt.
Kolom 1 heeft een negatieve coëfficiënt van de doelfunctie van -300, en komt als pivotkolom k in aanmerking.
Berekening van de coëfficiënten van deze kolom 1:
rij 1: 170 / 1 = 170 rij 2: 150 / 1 = 150 rij 3: a31 = 0, daar geen quotiënt berekenbaar.
Het kleinste quotiënt 150 krijgt men dus in rij 2. Met pivotelement a2,1 berekent men de nieuwe waarde van de doelfunctie als
- .
Ook kolom 2 heeft met -500 een negatieve coëfficiënt in de doelfunctie, en komt eveneens als pivotkolom k in aanmerking.
De berekeningen van de coëfficiënten van deze kolom 2:
rij 1: 170 / 2 = 85 rij 2: 150 / 1 = 150 rij 3: 180 / 3 = 60
Het kleinste quotiënt 60 wordt gevonden in rij 3. Met het mogelijke pivotelement a3,2 berekent men de nieuwe waarde van de doelfunctie als
- .
Het is dus gunstiger om als pivotelement het eerste getal uit rij r = 2 en kolom k = 1 te kiezen, aangezien men daarmee de doelfunctie het meest doet toenemen, namelijk van -36000 naar 9500.
Uitvoeren van een verwisselingsstap
[bewerken | brontekst bewerken]In de verwisselingsstap wordt de basisvariabele YB met de niet-basisvariabele x1 verwisseld.
Het simplex-tableau wordt volgens de bovenstaande regels omgerekend.
Het pivotelement is a2,1 = 1.
Berekening van de pivotelementen: a2,1 = 1 / 1 = 1 (met formule (1))
Berekening van de pivotrij r = 2: Elk element wordt door het pivotelement gedeeld. Kolom 1 bevat het pivotelement, en werd zonet berekend. Kolom 2: a2,2 = 1 / 1 = 1 (met formule (2)) rechterlid: b2 = 150 / 1 = 150 (met formule (3))
Berekening van de pivotkolom k = 1: Elk element wordt door het pivotelement gedeeld, het teken wordt omgekeerd. Doelfunctie: - (-300 / 1) = 300 Rij 1: a1,1 = - (1 / 1) = -1 (met formule (4)) Rij 2 is de pivotrij en werd reeds berekend Rij 3: a3,1 = (- (0 / 1) = 0
Berekening van de overige waarden - kolom 2, aangezien kolom 1 reeds berekend is: Doelfunctie kolom 2: -500 - (-300 × 1) /1 = -200 Rij 1 kolom 2: a1,2 = 2 - 1×1 / 1 = 1 (met formule (5)) Rij 2 kolom 2 behoort tot de pivotrij en is reeds berekend. Rij 3 kolom 2: a3,2 = 3 - 0×1 / 1 = 3
rechterlid: doelfunctie: G = -36000 - (-300 × 150) / 1 = 9500
Rij 1: b1 = 170 - 1×150 / 1 = 20 (met formule (6)) Rij 2 behoort tot de pivotrij en is al berekend Rij 3: b3 = 180 - 0×150 / 1 = 180
Het nieuwe simplex-tableau
[bewerken | brontekst bewerken]Na omrekening verkrijgt men een nieuw simplex-tableau:
-------------- ------------ | YB x2 | rechterlid ---- -------------- ------------ G | 300 -200 | 9500 ---- -------------- ------------ YA | -1 1 | 20 = b1 x1 | 1 1 | 150 = b2 YC | 0 3 | 180 = b3
De variabele x1 is tot de basis toegetreden, de variabele YB verschijnt in de titelrij en is dus een niet-basisvariabele.
Deze oplossing betekent: er worden 150 eenheden van product 1 gemaakt (rechterlid van de rij met x1). Van product 2 wordt niets vervaardigd (x2 is een niet-basisvariabele). Daarmee heeft het bedrijf een opbrengst van 9500 euro. Machine A staat 20 uur per maand stil (draait slechts 150 van de 170 mogelijke uren). Machine B draait op volle toeren. Machine C staat 180 uur stil, ze wordt dus helemaal niet gebruikt. Dit kan men ook afleiden uit de opgave: machine C wordt enkel bij de productie van product 2 gebruikt. Aangezien dit product niet vervaardigd wordt, gebruikt men machine C dan ook niet.
Verdere verbetering van de oplossing
[bewerken | brontekst bewerken]Aangezien de doelfunctie in het nieuwe simplex-tableau nog een negatieve coëfficiënt heeft, kan men die oplossing verder verbeteren. Dit gebeurt door een volgende simplex-iteratie.
Bij de keuze van de pivotelementen komt enkel kolom k = 2 in aanmerking, aangezien enkel hier de coëfficiënt van de doelfunctie negatief is. De pivotrij kiest men volgens het minimum van de coëfficiënten "rechterlid/kolom 2", dus het minimum uit
20 / 1 = 20 150 / 1 = 150 180 / 3 = 60
Men kiest dus rij r = 1 als pivotrij. Het pivotelement is a1,2 = 1, en de basisvariabele YA wordt met de niet-basisvariabele x2 verwisseld.
Na analoge berekeningen als hierboven wordt het nieuwe simplex-tableau:
-------------- ------------ | YB YA | rechterlid ---- -------------- ------------ G | 100 200 | 13000 ---- -------------- ------------ x2 | -1 1 | 20 x1 | 2 -1 | 130 YC | 3 -3 | 120
De optimale oplossing
[bewerken | brontekst bewerken]Aangezien de doelfunctie geen negatieve coëfficiënten meer bevat, is de optimale oplossing bereikt.
Er worden 130 eenheden van product 1 en 20 eenheden van product 2 gefabriceerd. Daarmee behaalt het bedrijf een winst van 13.000 euro. Machine A en B worden ten volle benut. Machine C draait 60 uur, en heeft dus nog een onbenutte capaciteit van 120 uur (in het tableau is YC = 120).
Opmerking over de keuze van het pivotelement
[bewerken | brontekst bewerken]In het voorbeeld wordt het pivotelement zodanig gekozen, dat de waarde van de doelfunctie bij de iteratiestap zo groot mogelijk wordt. Dit noemt men de Greatest-Change-Methode.
Alternatief kan men ook die kolom als pivotkolom nemen, waarvoor de coëfficiënt van de doelfunctie negatief en in absolute waarde het grootst is. Aangezien deze coëfficiënt een soort toename voorstelt, spreekt men van de Steepest-Unit-Ascent-Methode.
Vaak - maar niet altijd - levert de Greatest-Change-Methode met weinig simplex-iteraties een optimale oplossing.
Grafische oplossing
[bewerken | brontekst bewerken]In dit eenvoudige voorbeeld kent het probleem twee variabele parameters, in dat geval kan de oplossing grafisch voorgesteld worden in een 2-dimensionale figuur. Op de figuur worden de ongelijkheden voorgesteld als halfvlakken. De grensrechte wordt telkens afgebeeld. De doorsnede van deze halfvlakken bepaalt de regio met geldige oplossingen. Dit gebied is op de figuur geel gekleurd. De doelfunctie vormt ook steeds een rechte waarop het objectief eenzelfde waarde heeft, op de figuur zijn de rechten weergeven voor G = 0 en G = 13.000, in feite kan men een continue verzameling rechten tekenen. Bemerk dat al deze rechten die de doelfunctie voorstellen evenwijdig zijn, hoe groter G, hoe verder de rechte van de oorsprong ligt in dit geval. Men kan dan ook besluiten dat men zal wil nastreven deze rechte zo ver mogelijk van de oorsprong te krijgen, maar de rechte moet wel nog een punt binnen de geldige regio hebben, aangezien aan de nevenvoorwaarden moet voldaan zijn. Aangezien deze regio een convexe figuur is, ligt het maximum op een van de hoekpunten, dit is het snijpunt met de rechte waar G = 13.000, waar men x1 = 130 en x2 is 20 kan aflezen.
Op de figuur valt ook de werking van de simplex-iteraties grafisch te volgen. Het algoritme begon met beginoplossing (x1=0, x2=0), wat overeenkomt met het punt A. Dit is tevens het snijpunt van de rechten bepaald door de niet-negativiteitsvoorwaarden, namelijk de rechten met vergelijkingen x1=0 en x2=0. In tweede iteratie werd x1 aan de basis toegevoegd, YB werd nu een niet-basisvariabele. Wanneer men nu het snijpunt neemt van de rechten behorend bij deze niet-basisvariabelen, namelijk de rechte x2 = 0 en de rechte B, dan bekomt men het punt B(150,0), en dit blijkt inderdaad ook de oplossing te zijn die men vindt bij de tweede iteratiestap. Bemerk dat vanuit punt B, men de doelfunctie enkel kan doen toenemen wanneer men zich verplaatst langs rechte B, wanneer men zich verplaatst langs de x-as, daalt de doelfunctie. Het is dan ook inderdaad YB die uit de basis gehaald werd. In de derde iteratiestap zijn dan ook YA en YB niet-basisvariabelen: het snijpunt van de bijhorende rechten is opnieuw de oplossing van deze iteratiestap. Op deze manier loopt het simplex-algoritme langs de zijden van het geldig gebied, zolang er een zijde beschikbaar is langswaar de doelfunctie stijgend is, tot een maximum bereikt is.
Alternatief tableau
[bewerken | brontekst bewerken]Het simplex-tableau kan op een alternatieve manier ingedeeld worden, waarbij men alle variabelen (zowel basis- als niet-basisvariabelen) in de hoofding plaatst. Daaronder noteert men dan de doelfunctie en de nevenvoorwaarden, komt een variabele niet voor, dan schrijft men er een 0.
Het oorspronkelijke tableau ziet er dan als volgt uit:
----------------------------------- ------------- | x1 x2 YA YB YC | rechterlid ---- ----------------------------------- ------------- G | -300 -500 0 0 0 | ---- ----------------------------------- ------------- YA | 1 2 1 0 0 | 170 = b1 YB | 1 1 0 1 0 | 150 = b2 YC | 0 3 0 0 1 | 180 = b3
We zien dat de rijen/kolommen Ya,b,c een eenheidsmatrix vormen. Deze tabel toont de oplossing { x1 = 0, x2 = 0, YA = 170, YB = 150, YC = 180 } . Analoog aan voorheen kiest men kolom 2, omdat -500 in absolute waarde het grootst is, en een wijziging van x2 dus een potentieel grotere verandering in de objectieffunctie kan veroorzaken. Men kiest als pivotrij, rij YC. Dit besluit men opnieuw uit de berekening van rechterlid/x2. Voor rij YA bekomt men 170/2 = 85, voor rij YB 150/1 = 150, voor rij YC 180/3 = 60. Men ziet echter dat wanneer men voor x2 de waarde 85 zou gebruiken, de nevenvoorwaarde van rij YC nooit meer voldaan zou kunnen worden. Immers, vult men x = 85 in in 3 x2 YC = 255 YC, dan ziet men dat deze som nooit kleiner kan worden dan 180, aangezien YC positief is, en de ongelijkheid zou dus nooit voldaan kunnen worden. Daarom kiest men die rij waar de verhouding rechterlid/x2 het kleinst is, of dus rij YC. Ook hier voert men nu elementaire rijbewerkingen op de matrix uit, wat leidt tot het volgende tableau:
----------------------------------- ------------- | x1 x2 YA YB YC | rechterlid ---- ----------------------------------- ------------- G | -300 0 0 0 166,66 | ---- ----------------------------------- ------------- YA | 1 0 1 0 -0,666 | 50 YB | 1 0 0 1 -0,333 | 90 x2 | 0 1 0 0 0,333 | 60
Opnieuw is te zien dat x2 nu een basisvariabele is. De eenheidsmatrix is ook weer terug te vinden, met Ya,b en x2 In de voorlopige oplossing heeft x2 de waarde 60. De eerste kolom heeft nu een negatief element, rij YA heeft de kleinste kolomverhouding, pivoteren levert:
----------------------------------- ------------- | x1 x2 YA YB YC | rechterlid ---- ----------------------------------- ------------- G | 0 0 300 0 -33,33 | ---- ----------------------------------- ------------- x1 | 1 0 1 0 -0,666 | 50 YB | 0 0 -1 1 0,333 | 40 x2 | 0 1 0 0 0,333 | 60
Aangezien er nog één kolom met negatieve coëfficiënt blijft, levert een laatste maal pivoteren met rij YB:
----------------------------------- ------------- | x1 x2 YA YB YC | rechterlid ---- ----------------------------------- ------------- G | 0 0 200 100 0 | ---- ----------------------------------- ------------- x1 | 1 0 -1 2 0 | 130 YC | 0 0 -3 3 1 | 120 x2 | 0 1 1 -1 0 | 20
Men leest uiteindelijk dezelfde oplossing x1 = 130 en x2 = 20 af.
Algoritme in pseudocode
[bewerken | brontekst bewerken]In deze code wordt de begintabel als de matrix "tab" voorgesteld, waarbij G de onderste rij vormt en niet de bovenste. De matrix ziet er voor het voorbeeld als volgt uit:
----- ----- ------ | 1 | 2 | 170 | ----- ----- ------ | 1 | 1 | 150 | ----- ----- ------ | 0 | 3 | 180 | ----- ----- ------ |-300 |-500 |-36000| ----- ----- ------
Programma simplex() int r, k, i, j, pivot Invoer: "rij?", r Invoer: "kolom?", k pivot = tab[r,k] Voor i van 1 tot aantalRijen(tab) Voor j van 1 tot aantalKolommen(tab) Als i != r en j != k dan tab[i,j] = tab[i,j]-tab[i,k]*tab[r,j]/pivot EindeAls EindeVoor EindeVoor Voor j van 1 tot aantalKolommen(tab) tab[r,j] = tab[r,j]/pivot EindeVoor Voor i van 1 tot aantalRijen(tab) Als i != r dan tab[i,k] = 0 EindeAls EindeVoor Uitvoer: tab EindeProgramma