Voorraadtheorie. Mathijs van der Vlies 28 augustus Begeleiding: prof.dr. R. Núñez Queija

Voorraadtheorie Mathijs van der Vlies 28 augustus 2014 Begeleiding: prof.dr. R. N´ un ˜ez Queija

Korteweg-De Vries Instituut voor Wiskunde Faculteit der Natuurwetenschappen, Wiskunde en Informatica Universiteit van Amsterdam

Samenvatting Deze scriptie behandelt stochastische voorraadmodellen, waarbij kosten van een bedrijf moeten worden geminimaliseerd onder een onzekere vraag door middel van het goed kiezen van bestelmomenten en bestelgroottes. We behandelen eerst vastgestelde voorraadmodellen zoals Economic Order Quantity en het (s, Q)-naleveringsmodel, met bekende optimalisatiemethoden. Dan wordt een eigen model geconstrueerd, met exponentieel verdeelde levertijden en de restrictie dat er één bestelling tegelijk geplaatst kan worden. Er is dan sprake van een continue-tijds Markovketen. Dit model wordt eerst voornamelijk analytisch behandeld als (s, Q)-naleveringsmodel door middel van de matrix-geometrische methode. Hierbij wordt voor verschillende combinaties van het bestelpunt en de bestelgrootte de kosten op de lange termijn in kaart gebracht. De optimalisatie vindt plaats via een numerieke methode. Het model wordt behandeld als Markov beslissingsproces, waarna het algoritme van Successieve Approximatie toegepast wordt om de optimale bestelstrategie te vinden.

Titel: Voorraadtheorie Auteur: Mathijs van der Vlies, [email protected], 10251626 Begeleiding: prof.dr. R. N´ un ˜ez Queija Tweede beoordelaar: dr. N. Walton Einddatum: 28 augustus 2014 Korteweg-De Vries Instituut voor Wiskunde Universiteit van Amsterdam Science Park 904, 1098 XH Amsterdam http://www.science.uva.nl/math

2

Inhoudsopgave 1. Inleiding

4

2. Elementaire modellen 2.1. Economic Order Quantity . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2. Het (s, Q)-naleveringsmodel . . . . . . . . . . . . . . . . . . . . . . . . .

6 6 8

3. Het 3.1. 3.2. 3.3. 3.4. 3.5.

(s, Q)-model als Markovketen Formulering Markovproces . . . . . . . Matrix-geometrische methode . . . . . Voorraadkosten en kosten bij tekorten . Vaste kosten . . . . . . . . . . . . . . . Combinaties van s en Q . . . . . . . .

. . . . .

4. Markov beslissingsproces voor voorraad 4.1. Markov beslissingstheorie . . . . . . . . . 4.2. Successieve approximatie . . . . . . . . . 4.3. Formulatie MDP . . . . . . . . . . . . . 4.4. Truncatie van de toestandsruimte . . . . 4.5. Optimale bestelstrategie middels SA . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

17 18 20 27 28 29

. . . . .

32 32 34 39 42 43

5. Discussie

46

Bibliografie

47

A. SA-algoritme in MATLAB

48

B. Populaire samenvatting

51

3

1. Inleiding Voorraadmanagement is een belangrijk onderdeel van wat een bedrijf bezighoudt. Elk bedrijf dat producten verkoopt, zij het aan klanten of aan andere bedrijven, heeft ermee te maken. Voor elk product dat een bedrijf verkoopt, is er een voorraad van dit product. Slecht beheerde voorraad leidt tot tekorten, met als gevolg klantverlies. Het is voor veel bedrijven van uiterst belang dat de voorraad niet opraakt, omwille van service aan de klant. Goede service leidt tot een goede reputatie, wat weer leidt tot winst. Maar bedrijven kunnen hier ook in doorslaan. Het houden van zeer grote voorraden draagt ook weer kosten met zich mee. Denk aan verzekeringen voor de in voorraad gebrachte producten. Deze kosten vormen een significant deel van de kostprijs van de producten, significant genoeg om de voorraad zo laag mogelijk te houden. De crux van voorraadmanagement is dat de vraag naar het product vrijwel altijd onzeker is. Tekorten kunnen zonder waarschuwing ontstaan. Met dit in gedachten, heeft een bedrijf twee belangrijke vragen voor het houden van de voorraad: • Wanneer koopt het bedrijf opnieuw in? • Hoeveel koopt het bedrijf in? Er is geen snel antwoord te geven op deze vragen. Te grote bestelhoeveelheden leiden tot hoge voorraadkosten. Echter, wanneer de bestelhoeveelheden te klein zijn, moet het bedrijf te vaak bestellen, wat leidt tot oplopende vaste kosten (denk aan benzinekosten, administratiekosten, setupkosten bij productie). Kiest het bedrijf ervoor laat in te kopen, ontstaat het risico op tekorten tijdens de levertijd. Te vroeg inkopen leidt weer tot hoge voorraadkosten. Het oplossen van de bovenstaande vragen is waar de voorraadtheorie zich mee bezighoudt. Deze tak van wiskunde behoort tot de stochastische optimalisatie. Stochastische processen zoals Markovprocessen spelen hierbij een grote rol. De voorraadtheorie heeft als vakgebied veel raakvlakken met wachtrijtheorie, waarin ook regeneratieve processen voorkomen (processen die zich gedragen in cycli). Voorraadtheorie is een zeer toegepaste tak van wiskunde. Ze wordt altijd beoefend aan de hand van een model dat zo goed mogelijk de werkelijkheid probeert na te bootsen. Er bestaan al vele bekende modellen binnen de voorraadtheorie. De meeste van deze modellen staan erom bekend dat deze wiskundig zeer goed beheersbaar zijn. Voorbeelden hiervan zijn het EOQ-model en het (s, Q)-naleveringsmodel. In bedrijven worden deze vaak gebruikt, omdat ze resultaten leveren die vrijwel direct inzetbaar zijn. De simpliciteit van deze modellen heeft echter tot gevolg dat de gevonden optima niet altijd even betrouwbaar zijn. Het doel van dit onderzoek is een model te construeren dat relatief praktisch is en toch wiskundig exact kan worden opgelost.

4

De opbouw van de scriptie is als volgt. In het eerste hoofdstuk worden bekende modellen behandeld zoals de bovengenoemde. Het EOQ-model is hierbij een deterministisch voorraadmodel waarvoor de exacte oplossing eenvoudig te vinden is. Het (s, Q)naleveringsmodel wordt opgelost middels een heuristische methode, maar deze geeft geen exacte oplossing. Het tweede hoofdstuk is gewijd aan het aanpassen van de aannames van het (s, Q)naleveringsmodel, zodat het voorraadproces een Markovketen vormt. Zo maken we het model rijp voor een exacte analyse. De zogenaamde matrix-geometrische methode stelt ons in staat dit probleem exact te behandelen. Het laatste hoofdstuk laat de restrictie los dat elke bestelling dezelfde bestelgrootte moet hebben. We behandelen het resulterende model middels Markov beslissingstheorie. Deze theorie levert krachtige algoritmes voor het numeriek vinden van een optimale bestelstrategie. Er wordt een MATLAB-functie geconstrueert dat een van de algoritmes implementeert en zodoende de optimale bestelstrategie vindt. De lezer wordt aangeraden stof met betrekking tot Markovketens en lineaire algebra op elementair niveau goed te beheersen. Deze vakgebieden zullen veelvuldig worden gebruikt. Mijn dank gaat uit naar Sindo N´ un ˜ez Queija, voor het waardevolle advies dat hij heeft gegeven tijdens onze afspraken, voor de tijd die hij genomen heeft voor mij om dit project tot een geslaagd einde te brengen, en voor het voorstellen van een onderwerp in een vakgebied dat mij zal blijven interesseren.

5

2. Elementaire modellen Voor het opbouwen van de voorraadtheorie voor praktische doeleinden, zullen we verschillende voorraadmodellen beschouwen. Hierbij maken we ten eerste het onderscheid tussen deterministische en stochastische modellen. Het aspect van het model dat ofwel deterministisch, ofwel stochastisch wordt beschouwd, is de vraag naar het product. In het deterministische model is deze vraag bekend. In het stochastische model is deze niet a priori bekend, maar wordt de vraag gezien als een stochastische variabele, met een bekende kansverdeling. De nadruk zal in dit onderzoek liggen op stochastische modellen, maar het is toch nuttig om een deterministisch model te bekijken, omdat dit resultaten oplevert die ook voor stochastische modellen nuttig zijn. We bekijken een deterministisch voorraadmodel dat nuttige resultaten oplevert, het EOQ-model.

2.1. Economic Order Quantity Het EOQ-model is het simpelste model dat er is in de voorraadtheorie. Het model heeft een aantal aannames: 1. De vraag is deterministisch en constant in de tijd en we beschouwen de voorraad als een continue variabele 2. Aan de vraag moet direct voldaan worden 3. De aanvulorder mag willekeurig groot worden 4. De benodigde tijd om de voorraad aan te vullen is verwaarloosbaar 5. De gehele aanvulorder wordt op hetzelfde tijdstip geleverd 6. De enige relevante kosten zijn bestel- en voorraadkosten Uiteraard zijn veel van deze aannames in de praktijk niet realistisch, maar de resultaten van dit model zijn toch bruikbaar voor ingewikkeldere modellen. Herinner je dat de voorraadtheorie de volgende vragen beantwoordt: 1. Wanneer bestellen? 2. Hoeveel bestellen?

6

De eerste vraag is in dit model gemakkelijk te beantwoorden. Aangezien de vraag deterministisch is, er direct aan de vraag voldaan moet worden en de benodigde tijd om de voorraad aan te vullen verwaarloosbaar is, kunnen we simpelweg stellen dat we moeten bestellen wanneer de voorraad op is. De enige relevante vraag in dit model is dus wat de optimale bestelgrootte Q is. Hiertoe minimaliseren we de totale kosten T K(Q). Zoals gebruikelijk in voorraadmodellen, bekijken we alle kosten op jaarbasis. We bekijken dus de totale jaarlijkse kosten T K(Q). Deze valt uiteen in de volgende kosten: • K = de vaste kosten per order; • r = de jaarlijkse voorraadkosten per ge¨ınvesteerde euro in voorraad; • v = de inkoopprijs per eenheid product. Deze drie soorten kosten behoeven enige uitleg. Wat de vaste kosten K betreft, kun je denken aan transport- en administratiekosten. De jaarlijkse voorraadkosten moeten worden gezien als een percentage van het kapitaal dat in de voorraad wordt ge¨ınvesteerd. Dit zijn investeringen die niet tot opbrengsten leiden, in tegenstelling tot de inkoopkosten. Dit is geld dat bijvoorbeeld op de bank gezet had kunnen worden (een van de meest zekere vormen van winst). De winst die hierbij wordt misgelopen vormt over het algemeen het grootste deel van de voorraadkosten. Andere relevante kosten van deze soort zijn kosten van opslag en verzekering van de goederen. De resulterende kosten voor het in voorraad houden van één eenheid product gedurende één jaar zijn dus gelijk aan vr. Wat dit model karakteriseert is de deterministische vraag. In het EOQ-model is deze constant in de tijd; we geven de jaarlijkse vraag aan met D. Het aantal aanvulorders dat . Hieruit kunnen we de bestel- en inkoopkosten opmaken: dan gedaan wordt per jaar is D Q de jaarlijkse bestel- plus inkoopkosten =

D K + Dv. Q

Aangezien de vraag constant is, zal het voorraadniveau lineair dalen naar 0 en dan vanaf Q weer lineair dalen, ad infinitum. Het gemiddelde voorraadniveau is dus Q2 . Er geldt dus dat Q de jaarlijkse voorraadkosten = vr. 2 Dit geeft ons de volgende formule voor de totale jaarlijkse kosten T K(Q) =

D Q K + Dv + vr. Q 2

Het doel is nu deze kosten te minimaliseren. Dit is in dit geval gemakkelijk te doen door de afgeleide gelijk te stellen aan nul. Voor de optimale waarde Q∗ moet dus gelden: T K 0 (Q∗ ) = −

D 1 K + vr, 2 2 Q∗

en de oplossing r ∗

Q =

7

2KD vr

van deze vergelijking wordt de EOQ-formule genoemd. Om nog na te gaan dat dit een minimum is, nemen we de tweede afgeleide T K 00 (Q): T K 00 (Q) =

D K, 2Q3

en aangezien de te bestellen hoeveelheid Q altijd positief is, is de tweede afgeleide altijd positief. Dit bewijst dat de EOQ-formule inderdaad de jaarlijkse kosten minimaliseert. De EOQ-formule is een klassiek resultaat in de voorraadtheorie. Veel praktische situaties zijn met zekere aanpassingen terug te voeren op het EOQ-model, waardoor de EOQ-formule bruikbaar wordt. We zullen in het vervolg een bestelgrootte Q∗ die geoptimaliseerd is middels de EOQ-formule de economische bestelgrootte noemen. In praktische situaties zal de vraag naar een product niet met zekerheid bepaald kunnen worden; de vraag is stochastisch. Dit brengt extra problemen met zich mee. Bij een deterministisch voorraadmodel kon van te voren worden bepaald wanneer de voorraad op zou zijn, zodat nieuwe bestellingen precies aankwamen wanneer de voorraad op was. Daarentegen heb je bij een stochastisch model te maken met een kans op tekorten en overschotten. Bij tekorten kan een deel van de vraag niet direct uit voorraad geleverd worden; bij overschotten is er meer voorraad aanwezig dan vraag en worden er onnodige voorraadkosten gemaakt. Er zijn verschillende mogelijke redenen dat bedrijven tekorten willen beperken: • Er moet een bepaald serviceniveau behaald worden; het overgrote deel van de klanten moet direct uit voorraad zijn product geleverd krijgen. • Er is sprake van bederfelijke of seizoensgebonden goederen; producten die niet direct geleverd kunnen worden zijn waardeloos of kunnen slechts tegen een lagere prijs verkocht worden. • Wanneer een product niet direct geleverd kan worden, zijn klanten meer geneigd naar een concurrent over te stappen. Het is voor bedrijven ook gewenst overschotten te beperken: • De voorraadkosten zijn hoger wanneer een nieuwe bestelling aankomt terwijl de voorraad nog niet op is. • Er kan hier ook sprake zijn van bijvoorbeeld seizoensgebonden goederen; goederen die aan het einde van het seizoen nog niet verkocht zijn, kunnen slechts tegen een lagere prijs verkocht worden We behandelen nu een model dat een vergelijkbare opzet heeft met het EOQ-model. Het verschil is nu dat de vraag stochastisch is.

2.2. Het (s, Q)-naleveringsmodel De voorraadmodellen worden nu ingewikkelder. We hebben niet meer te maken met een deterministische vraag. In dit model wordt de voorraad beheerd over een lange

8

tijdsperiode. De parameters (kosten en verkoopprijzen) veranderen niet gedurende deze periode. De interessante toevoeging aan dit model is de stochastische vraag en positieve levertijden. Bij een deterministische vraag of afwezigheid van levertijd kon een bedrijf de tijd om te bestellen zo kiezen dat de bestelling aankwam zodra de voorraad op was. In dit model kan niet meer worden gegarandeerd dat een bestelling aankomt voordat de voorraad op raakt. Dit zorgt voor een extra beslissing in het model. Dit model behandelt dus twee keuzes: 1. Wanneer bestellen (bestelpunt s)? 2. Hoeveel bestellen (bestelgrootte Q)? Het bestelpunt wordt niet in tijdseenheden uitgedrukt, maar in eenheden product. Als de voorraad onder het bestelpunt is gedaald, vullen we de voorraad aan. Het bestelpunt zal doorgaans hoger liggen dan de verwachte vraag tijdens de levertijd. Immers, om goede service te bieden aan de klant moet er rekening gehouden worden met stochastische fluctuaties in de vraag. Een hoger bestelpunt kan worden gezien als een verzekering tegen het opraken van de voorraad tijdens de levertijd. De hoeveelheid product dat het bestelpunt boven de verwachte vraag tijdens de levertijd ligt, wordt de veiligheidsvoorraad genoemd. Een verhoging van de veiligheidsvoorraad verlaagt de kans dat het product uitverkocht raakt, maar zorgt voor een hoger gemiddeld voorraadsniveau. Een bedrijf moet een balans vinden tussen service aan de klanten en voorraadkosten. We kunnen twee mogelijke aannames doen over het optreden van vraag bij tekorten in een (s, Q)-model. 1. Nalevering. De vraag die optreedt als er geen voorraad aanwezig is, wordt geleverd zodra er weer voldoende voorraad is. 2. Verloren vraag. De vraag die optreedt als er geen voorraad aanwezig is, gaat verloren. We zullen in deze scriptie enkel het naleveringsmodel behandelen. We voeren de volgende terminologie in: 1. Voorraad op de planken. Dit is de voorraad die fysiek aanwezig is op de planken. Tekorten in voorraad worden hierbij niet meegerekend, dus deze voorraad is altijd niet-negatief. 2. Netto voorraad = (voorraad op de planken) - (na te leveren orders). Deze grootheid is alleen negatief als er nog naleveringen moeten plaatsvinden. In het model met verloren vraag is deze gelijk aan de voorraad op de planken. 3. Economische voorraad = (netto voorraad) + (orders in bestelling). Deze grootheid voegt aan de netto voorraad de nog niet binnengekomen aanvulorders toe. We zullen de voorraad beheren aan de hand van de economische voorraad. Als we alleen de netto voorraad zouden gebruiken, dan zouden we een aanvulorder kunnen doen terwijl er al een grote aanvulorder aankomt. Het model heeft de volgende aannames over het voorraadbeheer:

9

1. De voorraadpositie wordt continu bijgehouden en een aanvulorder kan op elk gewenst moment worden geplaatst. 2. De individuele vraagtransacties zijn zo klein dat het voorraadniveau gezien kan worden als een continue variabele. 3. Een aanvulorder ter grootte Q wordt geplaatst telkens als de economische voorraad tot het bestelpunt s daalt. 4. De levertijd van een bestelling is een positieve constante L∗ . 5. De gevraagde hoeveelheden in disjuncte tijdsintervallen zijn stochastisch onafhankelijk. Deze aannames zijn niet zo sterk als de aannames van de eerdere behandelde modellen. Dit model is dan ook in veel praktische situaties bruikbaar. We zullen het (s, Q) voorraadmodel niet exact behandelen, maar een aanpak geven die de exacte oplossing benadert. Een voordeel van deze aanpak is dat de stochastische vraag niet geheel beschreven hoeft te worden. Alleen de vraag gedurende de levertijd wordt bekeken: XL = de totale vraag gedurende de levertijd. Evenzo bekijken we de volgende grootheden: fL (x) = kansdichtheid van de vraag gedurende de levertijd µL = de verwachte waarde van de vraag gedurende de levertijd σL = de standaardafwijking van de vraag gedurende de levertijd In feite hebben we alleen µL en σL nodig voor de analyse van dit probleem. Deze worden geschat uit data over de vraag. Zoals hierboven al vermeld is, zullen we een heuristische analyse geven van het naleveringsmodel. De exacte analyse is nogal complex en we hebben hier te maken met een praktische tak van wiskunde. In de praktijk hebben we slechts beperkte informatie (data) over het stochastische gedrag van de vraag naar een product. Het is dan ook wenselijk dat de oplossing van het model slechts beperkte informatie nodig heeft. De heuristische oplossing voldoet hieraan. Voor een gegeven (s, Q) met s > 0 leiden we benaderingen af voor de volgende prestatiematen: • de gemiddelde voorraad op de planken • de gemiddelde achterstand in levering • de kans dat het systeem buiten voorraad raakt gedurende de levertijd • de fractie vraag die direct uit voorraad geleverd wordt.

10

Als we eenmaal uitdrukkingen hebben gevonden voor deze grootheden, kunnen we gemakkelijk uitdrukkingen formuleren voor de kosten die we vervolgens optimaliseren. We bepalen deze langetermijngemiddelden uit het gedrag van het voorraadsysteem gedurende één cyclus: het tijdsinterval tussen twee opeenvolgende tijdstippen waarop een aanvulorder binnenkomt. Aan de hand van dit begrip definiëren we vervolgens: I1 = E[voorraad op de planken aan het einde van een cyclus] I2 = E[voorraad op de planken aan het begin van een cyclus] en S1 = E[tekort aan het einde van een cyclus] S2 = E[tekort aan het begin van een cyclus]. Het volgende resultaat staat centraal in de analyse: de netto voorraad vlak voordat een aanvulorder binnenkomt = s − XL .

(2.1)

Dit resultaat is intu¨ıtief gemakkelijk te bewijzen. Allereerst merken we op dat de levertijd in dit model constant is. Dit betekent dat orders in dezelfde volgorde aankomen als dat ze worden besteld. De netto voorraad is gedefinieerd door de voorraad op de planken minus het aantal na te leveren orders. Op het moment dat een aanvulorder S geplaatst wordt, is de economische voorraad precies gelijk aan s. Aangezien de volgorde van bestellingen behouden blijft, zullen alle orders, die voor S geplaatst zijn, binnen zijn vlak voordat de order S aangekomen is. Aangezien tekorten worden nageleverd, volgt dat de netto vooraad vlak voordat S aangekomen is gelijk is aan s minus de totale vraag gedurende de levertijd, wat het bewijs voltooit. Deze uitdrukking stelt ons in staat de bovengenoemde grootheden uit te rekenen. Ten eerste worden I1 en S1 gegeven door I1 = E[s − XL ]+ ,

S1 = E[XL − s]+ .

Uit (2.1) leiden we direct af dat de netto voorraad vlak na binnenkomst van een aanvulorder gelijk is aan s + Q − XL . Dit geeft ons de volgende uitdrukkingen voor I2 en S2 : I2 = E[s + Q − XL ], S1 = E[XL − s − Q]+ . Deze verwachtingen kunnen we middels de Law of the Unconcious Statistician weer schrijven in termen van de kansdichtheid fL van XL : Z s Z s+Q I1 = (s − x)fL (x)dx, I2 = (s + Q − x)fL (x)dx 0

en

0

∞

Z

Z (x − s)fL (x)dx,

S1 =

0

0

(x − s − Q)fL (x)dx.

S2 =

s

De relatie Z a Z (a − x)fL (x)dx =

∞

s+Q

∞

Z (a − x)fL (x) −

∞

Z (a − x)fL (x) = a − µL +

a

(x − a)fL (x)dx, a

11

∞

geeft ons vervolgens het volgende verband tussen I1 en S1 respectievelijk I2 en S2 : I1 = s − µL + S1 ,

I2 = s + Q − µL + S2 .

(2.2)

De uitdrukkingen voor S1 en S2 kunnen weer verder vereenvoudigd worden wanneer we gaan kijken naar specifieke verdelingen voor XL , zoals we later in een voorbeeld terug zullen zien. Middels deze uitdrukkingen kunnen we de bovengenoemde prestatiematen benaderen. Gemiddelde voorraad op de planken Door niet te kijken naar de lange termijn voorraad, maar naar de voorraad gedurende een cyclus, verkrijgen we de volgende benadering voor de gemiddelde voorraad op de planken: 1 de gemiddelde voorraad op de planken ≈ (I1 + I2 ). 2 Deze uitdrukking kunnen we middels (2.2) herschrijven tot 1 1 de gemiddelde voorraad op de planken ≈ s − µL + Q + (S1 + S2 ) 2 2

(2.3)

Gemiddelde achterstand in levering Analoog aan het bovenstaande gemiddelde benaderen we deze door 1 de gemiddelde achterstand in levering ≈ (S1 + S2 ). 2

(2.4)

De kans op buiten voorraad raken Vanwege de formule voor de netto voorraad vlak voor binnenkomst van een bestelling (2.1) is het redelijk de kans dat het systeem buiten voorraad raakt gedurende de levertijd te benaderen door P (XL > s), dat wil zeggen: Z ∞ de kans op buiten voorraad raken gedurende de levertijd ≈ fL (x)dx. (2.5) s

Fractie direct geleverde vraag Deze grootheid is zeer belangrijk voor het bepalen van het bestelpunt s. De fractie direct geleverde vraag ten opzichte van de totale vraag moet gemaximaliseerd worden voor een goede service aan de klant. De fractie direct geleverde vraag kan worden uitgedrukt door de volgende stochasten. Definieer de stochasten D(t) en V (t) als de totale hoeveelheid vraag in het tijdsinterval (0, t) respectievelijk de totale hoeveelheid vraag in (0, t) waaraan niet direct voldaan kan worden. De fractie van de vraag, waaraan op de lange termijn niet direct voldaan kan

12

worden, is dan gelijk aan de limiet limt→∞ V (t)/D(t). Deze limiet hangt samen met de verwachtingen voor beide stochasten. Namelijk, de limiet is gelijk aan E[hoeveelheid vraag per cyclus die niet direct leverbaar is] . E[totale vraag in een cyclus]

(2.6)

Deze formule is intu¨ıtief duidelijk, maar het bewijs vergt geavanceerde kansrekening die in deze scriptie niet aan bod komt. Uit de definitie van S1 en S2 volgt dat de teller van (2.6) gelijk is aan E(S1 − S2 ). Ook de noemer is makkelijk te vinden. Merk hiertoe op dat in het naleveringsmodel uiteindelijk aan alle vraag is voldaan. Op de lange termijn zal de gemiddelde vraag per cyclus dus gelijk zijn aan de gemiddelde voorraad die besteld wordt per cyclus, wat gelijk is aan Q. Oftewel, E[totale vraag in een cyclus] = Q. Door nu de integraalrepresentaties voor S1 en S2 te gebruiken, vinden we het volgende belangrijke resultaat: de fractie vraag die niet direct uit voorraad geleverd wordt 1 = Q

Z

∞

Z

∞

(x − s − Q)fL (x)dx .

(x − s)fL (x)dx −

(2.7)

s+Q

s

We zullen nu deze prestatiematen gebruiken om een optimale waarde te vinden voor (s, Q). We maken hierbij een afweging tussen het minimaliseren van tekorten en het minimaliseren van voorraad- en bestelkosten. Omdat het vaak lastig is het effect van tekorten in de voorraad te kwantificeren, wordt er meestal gebruik gemaakt van een service-eis. We behandelen deze aanpak in de volgende sectie. Minimalisering van kosten onder een service-eis Net als in het EOQ-model nemen we aan dat de volgende gegevens van de kostenstructuur bekend zijn: K = vaste kosten verbonden aan een aanvulorder r = voorraadkosten per in voorraad ge¨ınvesteerde geldeenheid per tijdseenheid v = inkoopkosten per eenheid Deze kosten zullen worden geminimaliseerd onder een service-eis. Een service-eis neemt vaak een van de volgende twee vormen aan: P1 De kans op geen tekort gedurende de levertijd van een aanvulorder moet minstens gelijk zijn aan een zelfgekozen getal α met 0 < α < 1. P2 De fractie van de vraag die direct uit voorraad geleverd wordt moet ten minste gelijk zijn aan een gegeven getal β met 0 < β < 1.

13

We zullen in een voorbeeld zien dat de tweede servicemaat geschikter is om de service aan de klant te meten dan de eerste. Voor het minimaliseren van de kosten onder de service-eisen P1 of P2 kijken we naar de gemiddelde kosten per tijdseenheid. We hebben al de nodige uitdrukkingen afgeleid om deze kosten te bepalen. Namelijk, de gemiddelde voorraad op de planken wordt gegeven door (2.3), terwijl de gemiddelde bestelkosten per tijdseenheid volgen met de formule het gemiddelde aantal aanvulorders per tijdseenheid = µ1 /Q,

(2.8)

waarbij µ1 de gemiddelde vraag per tijdseenheid voorstelt. De gemiddelde kosten per tijdseenheid worden geminimaliseerd onder een van de restricties P1 of P2 . Wiskundig gezien zouden we nu s en Q simultaan moeten bepalen. De stochastische aard van de vraag maakt dit echter een bewerkelijk proces. Het blijkt dat in de praktijk een veel eenvoudigere aanpak vrijwel even goed werkt. Deze aanpak berust erin eerst de bestelgrootte Q te bepalen en daarna het bestelpunt s.

2.2.1. De sequenti¨ ele aanpak Eerst gebruiken we de EOQ-formule om de ’optimale’ bestelgrootte Q0 te bepalen: r 2µ1 K . (2.9) Q0 = vr Voor het bepalen van het optimale bestelpunt s bij deze bestelgrootte wenden we ons tot de prestatiematen waarvoor we uitdrukkingen hebben gegeven. Middels de uitdrukking (2.5) wordt de service-eis P1 gegeven door Z ∞ fL (x)dx ≤ 1 − α. s

De service-eis P2 wordt met behulp van (2.7) gegeven door Z ∞ Z ∞ 1 (x − s)fL (x)dx − (x − s − Q0 )fL (x)dx ≤ 1 − β. Q0 s+Q0 s Aangezien beide prestatiematen stijgend zijn in s, evenals de voorraad- en bestelkosten, worden deze service-eisen verzadigd bij een optimaal bestelpunt. Voor het bepalen van het bestelpunt lossen we dus de voor P1 en P2 een van de volgende vergelijkingen op: Z ∞ fL (x)dx = 1 − α, (2.10) s

respectievelijk 1 Q0

Z

∞

Z

∞

(x − s)fL (x)dx − s

(x − s − Q0 )fL (x)dx

= 1 − β.

(2.11)

s+Q0

De sequentiële aanpak blijkt uitstekend te werken wanneer de economische bestelgrootte Q0 voldoende groot is. Om voldoende rekening te houden met schommelingen in

14

de vraag hebben we voor deze aanpak nodig dat Q0 groter is dan de standaardafwijking van de vraag gedurende de levertijd. We hebben nu een vergelijking gevonden waaruit de optimale bestelgrootte afgeleid kan worden. Deze vergelijking is in het algemene geval nog zeer bewerkelijk: afhankelijk van de dichtheid fL (x) kan het oplossen van de vergelijking (2.11) nog zeer ingewikkeld zijn. In het geval van een normaal verdeelde vraag kunnen we de vergelijking echter verder vereenvoudigen. Optimaal bestelpunt voor normaal verdeelde vraag Uiteraard zal de verdeling van de vraag gedurende de levertijd in de praktijk niet precies bekend zijn. Vaak wordt dan de normale verdeling gebruikt om de vraag te modelleren. Dit wordt gerechtvaardigd als de vraag ontstaat uit een groot aantal afnemers vanwege de centrale limietstelling. Stel dus dat de vraag XL gedurende de levertijd normaal verdeeld is, met verwachting µL en standaardafwijking σL . Aangezien de vraag nooit negatief is, eisen we dat σL /µL niet te groot is (zeg, σL /µL ≤ 0.5). Voor het vereenvoudigen van de vergelijkingen (2.10) en (2.11) gebruiken we de volgende representatie voor het bestelpunt s: s = µL + kσL We zullen nu het getal k bepalen, waaruit we s afleiden. Dit getal wordt de veiligheidsfactor genoemd, omdat kσL (= s − µL ) de veiligheidsvoorraad weergeeft. Deze representatie voor s stelt ons in staat vergelijking (2.10) in termen van de standaardnormale verdeling te herschrijven. Aangezien X L − µL > k = 1 − Φ(k), P (XL > s) = P σL met Φ(k) de verdelingsfunctie van de standaardnormale verdeling, kunnen we de vergelijking (2.10) vereenvoudigen tot 1 − Φ(k) = 1 − α.

(2.12)

De factor k wordt dan berekend middels een numerieke benadering voor de inverse functie Φ−1 (k) Voor het oplossen van vergelijking (2.11) hebben we de normale verliesfunctie nodig. Deze functie wordt gedefinieerd door Z ∞ 1 2 1 I(z) = √ (x − z)e− 2 x dx. 2π z De functiewaarden worden numeriek bepaald, omdat de integraal niet exact op te lossen is. Om (2.11) te herschrijven middels I(z), merken we op dat de functie gelijk is aan de volgende verwachting: I(z) = E(U − z)+ ,

15

waarbij U een standaardnormaal verdeelde stochast is. Aangezien nu voor elk getal a, Z ∞ (x − a)fL (x)dx = E(XL − a)+ a " + # XL − µL a − µL = σL E − σL σL a − µL = σL I , σL kunnen we de integraalrepresentaties voor S1 en S2 vereenvoudigen tot s − µL s + Q − µL S1 = σ L I , S2 = σL I . σL σL

(2.13)

Dit geeft ons de volgende uitdrukking voor vergelijking (2.11): s + Q − µL s − µL + σL I = (1 − β)Q0 . σL I σL σL Om deze vergelijking nog verder te vereenvoudigen, merken we op dat de tweede term gedefinieerd is als het verwachte tekort vlak na de binnenkomst van een bestelling. Als het vereiste serviceniveau hoog genoeg is, zal het bijna nooit plaatsvinden dat er nog een tekort is vlak nadat een aanvulorder binnengekomen is. Voor praktische doeleinden wordt dan ook meestal de tweede term verwaarloosd. We hoeven dan alleen nog de simpele vergelijking σL I(k) = (1 − β)Q0 (2.14) op te lossen. Het getal k wordt dan berekend middels een numerieke benadering voor de inverse functie I −1 (k). Deze vergelijking is een van de meest belangrijke wetenschappelijke resultaten in de voorraadtheorie. Het geeft een eenvoudige manier om het optimale bestelpunt s te berekenen gegeven een bestelgrootte Q0 . Aangezien we hier met een zeer specifieke verdeling (de normale verdeling) werken, welke statistisch zeer goed hanteerbaar is, is deze vergelijking in de praktijk goed bruikbaar. Het geeft echter geen exacte oplossing voor een optimale combinatie van (s, Q). We hebben al enkele redenen genoemd waarom deze oplossing niet exact is, waarvan de belangrijkste zijn: • De grootheden s en Q worden niet simultaan berekend: het model bepaalt eerst Q en dan s. • De prestatiematen ’gemiddelde voorraad op de planken’, ’gemiddelde achterstand in levering’ en ’kans op buiten voorraad raken’ worden niet exact berekend. • De term S2 wordt verwaarloosd in de service-eis P2 . Het belangrijkste kritiekpunt is het niet simultaan berekenen van s en Q. Maar hoe kunnen we deze grootheden dan wel tegelijk optimaliseren? We zullen het gehele model herbouwen. We maken hierbij gebruik van continue Markovketens en brengen hiermee het verloop van de economische en fysieke/netto voorraad in kaart.

16

3. Het (s, Q)-model als Markovketen Om een exacte oplossing te geven voor (s, Q) moeten we een beeld hebben van het verloop van de voorraad gedurende de tijd. Markovketens zijn hiervoor uitermate geschikt. Voor het werken met Markovketens hebben we de volgende aannames nodig: 1. De voorraadpositie wordt continu bijgehouden en een aanvulorder kan op elk gewenst moment worden geplaatst. 2. De vraag verloopt volgens een Poissonproces: individuele vraagtransacties zijn precies gelijk aan 1 en de tijd tussen transacties is exp(λ) verdeeld. 3. Een aanvulorder ter grootte Q ∈ N wordt geplaatst telkens als de economische voorraad tot het bestelpunt s ∈ N daalt. 4. De levertijd van een bestelling is een stochastische variabele die exp( L1 ) verdeeld is. 5. De gevraagde hoeveelheden in disjuncte tijdsintervallen zijn stochastisch onafhankelijk (Markov-eigenschap). We leggen aan ons model verder nog een belangrijke restrictie op: er kan maar énén enkele bestelling tegelijk worden geplaatst. Deze situatie doet zich bijvoorbeeld voor wanneer een bedrijf één vrachtwagen tot zijn beschikking heeft dat op pad gestuurd wordt voor bestellingen. Deze aannames verschillen op een aantal punten van het normale (s, Q)-voorraadmodel. Een is dat de vraag niet meer een continue variabele is. De vraag verloopt nu in gelijke sprongen die onafhankelijk van elkaar zijn en waarbij de tijd tussen sprongen stochastisch is. Hiervoor hebben we een aftelbare toestandsruimte I = N voor de vraag nodig; de bestelgrootte Q en het bestelpunt s moeten dus ook natuurlijke getallen zijn. Een belangrijk tweede verschil is dat de levertijd niet meer constant is, maar een exponentieel verdeelde stochast. Deze aanname is zeer nuttig voor het uitwerken van het model. De exponentiële verdeling heeft de gunstige eigenschap dat deze geheugenloos is: wanneer de tijd dat een gebeurtenis plaatsvindt exponentieel verdeeld is, is het plaatsvinden van deze gebeurtenis onafhankelijk van de tijd dat er gewacht is. Voor de levertijd betekent dit dat het voor het bepalen van het leveringsmoment het niet van belang is te weten wanneer de bestelling geplaatst is, enkel dat deze geplaatst is. De kostenstructuur van dit model ziet er als volgt uit: • Voorraadkosten r per tijdseenheid per eenheid product in voorraad • Vaste kosten K per bestelling • Boetekosten b per tijdseenheid per eenheid product tekort

17

Er kan een probleem ontstaan in dit model. Het verbieden van bestellingen tijdens het verloop van een andere bestelling kan zorgen voor een niet-positief recurrent model. Namelijk als Q klein is en de verwachte vraag λ per tijdseenheid groot is, kunnen er tekorten in de voorraad ontstaan die niet meer nageleverd kunnen worden omdat de bestelde hoeveelheid per keer te klein is en er maar één bestelling tegelijk plaats kan vinden. We beschouwen daarom twee mogelijke strategieën voor het regelen van het voorraadniveau: • een (s, Q)-strategie: er wordt Q besteld wanneer de economische voorraad gedaald is tot s of lager en wanneer de vorige bestelling geleverd is; • een (s, S)-strategie: het verschil met een bepaald referentieniveau S voor de economische voorraad wordt hersteld wanneer de economische voorraad gedaald is tot s of lager. Aangezien er maar één bestelling tegelijk geplaatst kan worden, zal het voor kunnen komen dat de nettovoorraad zeer laag is op het moment dat een bestelling binnenkomt. Het kan dan wenselijk zijn extra te bestellen. We verwachten daarom dat de (s, S)strategie een superieure strategie zal zijn voor dit probleem. Toch zullen we voor inzicht in het voorraadproces in eerste instantie ons voornamelijk bezig houden met de bespreken van (s, Q)-strategieën. We zullen voor dit model het gedrag van de nettovoorraad in kaart brengen en zodoende uitdrukkingen geven voor de verscheidene vormen van kosten (vaste kosten, voorraadkosten, kosten bij tekorten). Wanneer we deze kosten hebben gevonden, kunnen we deze vervolgens in kaart brengen voor verschillende combinaties s en Q. Het feit dat de levertijd exponentieel verdeeld is, zal hiertoe een belangrijke rol spelen. Aangezien de vraag verder ook exponentieel verdeeld is, kunnen we het voorraadproces beschrijven middels een continue-tijds Markovketen (CTMC). Er bestaan vele technieken op deze processen die ons helpen exacte uitdrukkingen te geven voor de kosten op de lange termijn. We zullen eerst het betreffende CTMC formuleren voor ons probleem.

3.1. Formulering Markovproces Een continue-tijds Markovketen met aftelbare toestandsruimte kan worden gerepresenteerd door zijn genererende matrix Q. Deze is voor het model met exponentiëel verdeelde levertijden gemakkelijk te geven. We merken op dat er slechts één bestelling tegelijk kan worden geplaatst. Deze bestelling wordt geplaatst wanneer de nettovoorraad onder s ligt. We kunnen het voorraadproces opdelen in twee deelprocessen die beide CTMC’s zijn: het vraagproces en het bestelproces. Het vraagproces is simpelweg een sterfteproces: de nettovoorraad daalt met stappen van 1, waarbij de tijd tussen ’sterftes’ (de aankomst van klanten) exp(λ) verdeeld is. Analoog wordt het bestelproces gegeven door de nettovoorraad te laten stijgen met stappen van Q, waarbij de tijd tussen bestellingen exp( λ1 ) verdeeld is, met als voorwaarde dat de nettovoorraad onder s moet liggen. De aard van deze CTMC is zo dat er één bestelling tegelijk wordt geplaatst, analoog aan het feit dat klanten een voor een aankomen in het vraagproces. Deze twee processen tellen op tot

18

het voorraadproces, dat dus gegeven wordt door de genererende matrix Q = (qij )i,j≤s+Q , met qi,i−1 = λ ∀i, 1 i ≤ s, qi,i+Q = L ( −λ qii = −λ + L1 qij = 0

als i > s, als i ≤ s,

anders.

De eerste vraag die we ons bij dit proces moeten stellen is wanneer deze positief recurrent is. Voor dit specifieke proces is het zo dat, wanneer Q te klein is ten opzichte van λ, de nettovoorraad na lange tijd zal blijven dalen. Aangezien er maar één bestelling tegelijk geplaatst kan worden en de levertijd positief is, kan het zo zijn dat er na lange tijd meer klanten binnenkomen dan dat er van het product kan worden besteld. In deze situatie is er geen sprake van een evenwichtsverdeling, aangezien de nettovoorraad op de lange termijn naar −∞ gaat. We willen dit scenario voorkomen door een voldoende en noodzakelijke conditie vast te stellen zodat het proces positief recurrent is. Om deze conditie vast te stellen, geven we eerst de genererende matrix weer:   −λ λ 0 ··· 0 0 0 0 0 ···   .. .. .. .. ..   0 ... ... ... . . . . . · · ·   . .. .. .. .. .. .. ..  .  . 0 . . . . . . · · ·  .   . .. ..  . . λ . 0 0 0 0 · · · Q= . .   0 · · · · · · 0 −λ λ 0 0 0 · · ·   1   1 0 · · · · · · 0 −λ − λ 0 0 · · ·   L L 1 1   0 0 · · · 0 0 −λ − λ 0 · · ·   L L .. .. .. .. .. .. .. .. .. . . . . . . . . . . . . Deze matrix heeft een speciale vorm. We kunnen de matrix namelijk opdelen in vierkante matrices ter grootte Q, zodat deze als volgt gerepresenteerd kan worden:   B0 A0 0 0 0 · · · B1 A1 A0 0 0 · · ·    0 A2 A1 A0 0 · · ·   (3.1) Q =  0 0 A A A · · · , 2 1 0    0 0 0 A2 A1 · · ·   .. .. .. .. .. . . . . . . . . waarbij de nulentries vierkante nulmatrices ter grootte Q zijn en B0 , B1 , A0 , A1 en A2

19

matrices van dezelfde grootte, die gegeven worden door   −λ λ 0 ··· 0 ..   .. .. .. . . .  0 .   .  . . .  .. .. .. 0  B0 =  .. ,  .  ... ..  .. . λ 0 · · · · · · 0 −λ  1 0 ··· 0 L . . ..  . ..   0 .. B1 = A2 =  . . , .. ... 0  .. 0 · · · 0 L1   0 0 ··· 0  .. .. ..  . . . A0 =  ,  0 0 · · · 0 λ 0 ··· 0  0 ··· 0 −λ − L1 λ ..  .. .. .. . . .  0 .  . . . .  . . . . A1 =  . . . . 0  . ... ...  .. λ 0 · · · · · · 0 −λ −

(3.2)

(3.3)

(3.4)



1 L

   .   

(3.5)

Een matrix die op deze manier gepartitioneerd kan worden heet een matrix in matrixgeometrische vorm. Het blijkt dat oplossingen van de balansvergelijkingen van zulke matrices gegeven kunnen worden door het oplossen van een eindig stelsel vergelijkingen, in tegenstelling tot een oneindig stelsel waar we in eerste instantie mee te maken hebben. De methode waarop dit gebeurt, wordt beschreven in de volgende sectie.

3.2. Matrix-geometrische methode De vorm die de matrix Q aanneemt zegt veel over de structuur van het proces. We zien dit in wanneer we de toestanden i ≤ s + Q representeren als tweedimensionale toestandsvectoren (η, k), waarbij η ∈ H, en k ∈ K voor zekere aftelbare verzameling H en eindige verzameling K. Toestanden worden ingedeeld in niveaus η. De tridiagonale blokstructuur van de matrix zorgt er vervolgens voor dat overgangen tussen toestanden alleen als volgt plaats kunnen vinden: • tussen toestanden van hetzelfde niveau; • naar toestanden van één niveau hoger; • naar toestanden van één niveau lager.

20

Voor het behandelen van de matrix-geometrische vorm in het algemeen zullen we er in het vervolg voor het gemak van uit gaan dat H = {0, 1, . . . } en K = {1, 2, . . . , m} voor zekere m ∈ N. We zijn ge¨ınteresseerd in de evenwichtsverdeling van een dergelijke matrix in matrixgeometrische vorm. We groeperen hiertoe, analoog aan de toestanden, de evenwichtsverdeling in subvectoren ter lengte m door te nemen π = (π0 , π1 , . . . ), met πη = (π(η, 1), . . . , π(η, m)). Als we nu een matrix Q gebruiken van de vorm (3.1) samen met de gepartitioneerde πη , levert dit het volgende stelsel balansvergelijkingen op: π0 B0 + π1 B1 = 0 π0 A0 + π1 A1 + π2 A2 = 0 π1 A0 + π2 A1 + π3 A2 = 0, .. . πi−1 A0 + πi A1 + πi+1 A2 = 0 i = 2, 3, . . . Deze vergelijkingen zijn enkel op te lossen wanneer het proces positief recurrent is. De matrix-geometrische methode geeft een voldoende en noodzakelijke conditie voor positieve recurrentie. Stelling 3.1. Laat A = A1 + A2 + A3 en zij πA de evenwichtsverdeling bij A (deze bestaat aangezien alle rijen van A sommeren tot 0 en A een eindige-dimensionale matrix is). Dan is het irreducibele Markovproces positief recurrent dan en slechts dan als πA A0 e < πA A2 e,

(3.6)

waarbij e een vector van enen is met lengte m. Het bewijs van deze stelling maakt gebruik van genererende functies, die buiten het domein van dit onderzoek liggen. We verwijzen de ge¨ınteresseerde lezer naar Neuts [5]. De driftconditie heeft wel een intu¨ıtieve interpretatie. De linker- en rechterkant van de driftconditie zijn namelijk gelijk aan de gemiddelde overgangsintensiteiten van een niveau η ≥ 2 naar het niveau erboven en beneden respectievelijk. Volgens de driftconditie moet dus gelden dat de gemiddelde overgangsintensiteit naar boven kleiner is dan die naar beneden. Dit is precies wat nodig is voor positieve recurrentie (de niveaus zijn van boven onbegrensd en van beneden begrensd). Wanneer aan deze conditie voldaan is, kan de evenwichtsverdeling worden uitgerekend. Het fundamentele resultaat van de matrix-geometrische methode is nu dat er tussen opeenvolgende πη ’s een geometrisch verband bestaat, dat wil zeggen, πi = πi−1 R,

21

i = 2, 3, . . .

(3.7)

voor een zekere vaste vierkante matrix R ter grootte m. Inductief levert dit op: πi = π1 Ri−1 , i = 2, 3, . . .

(3.8)

Deze vergelijking impliceert dat we alle πi voor i = 2, 3, . . . kunnen bepalen middels R, π0 en π1 . Eerst behandelen we hoe R bepaald wordt. Deze vergelijking substitueren in het stelsel balansvergelijkingen levert ons op: π1 Ri−2 A0 + π1 Ri−1 A1 + π1 Ri A2 = 0, oftewel π1 Ri−2 (A0 + RA1 + R2 A2 ) = 0, wat leidt tot de karakteristieke vergelijking voor R: A0 + RA1 + R2 A2 = 0.

(3.9)

Deze vergelijking is vaak lastig exact op te lossen, maar er bestaat wel een numerieke methode die de juiste matrix R benadert. Ten eerste vermenigvuldigen we de vergelijking met A−1 1 om te krijgen: 2 −1 A0 A−1 1 + R + R A2 A1 = 0, waarbij we R naar de andere kant halen zodat 2 R = −A0 A1−1 − R2 A2 A−1 1 = −V − R W.

We benaderen R vervolgens met de volgende rij R(k) : R(0) = 0,

2 R(k) = −V − R(k−1) W,

k = 1, 2, . . .

(3.10)

Zo kunnen we R benaderen tot de gewenste precisie. Ten slotte dienen we nog π0 en π1 uit te rekenen. Hiertoe gebruiken we de eerste twee van de balansvergelijkingen. Deze luiden: π0 B0 + π1 B1 = 0, π0 A0 + π1 A1 + π2 A2 = 0. Door π2 te substitueren met π1 R kunnen we deze twee vergelijkingen in blokmatrixvorm schrijven: B0 A0 (π0 , π1 ) = (0, 0), (3.11) B1 A1 + RA2 waaruit we een oplossing π ˆ0 en π ˆ1 kunnen destilleren. Deze moeten nog genormaliseerd worden door de conditie πe = 1. Deze conditie bevat een oneindige som, maar is concreet te maken door gebruik te maken van de geometrische eigenschap: ∞ X 1 = πe = π0 e + πi e = π0 e + = π0 e +

i=1 ∞ X i=1 ∞ X

π1 Ri−1 e π1 R i e = π0 e + π1

i=0

∞ X i=0

22

! Ri

e

Om aan de conditie πe = 1 te voldoen moet de reeks bestaat de inverse (I − R)−1 en verkrijgen we

P∞

i=0

Ri convergeren. In dat geval

π0 e + π1 (I − R)−1 e = 1. Om dus aan π0 en π1 te komen, berekenen we α=π ˆ0 e + π ˆ1 (I − R)−1 e. De evenwichtskansen π0 en π1 worden dan gegeven door π0 =

π ˆ0 , α

π1 =

π ˆ1 . α

Daarmee hebben we gelijk de volledige evenwichtsverdeling door te gebruiken dat πi = π1 Ri−1 voor i = 2, 3, . . . .

3.2.1. Toepassing op het voorraadmodel We hadden al vastgesteld dat de genererende matrix van het voorraadprobleem de vorm (3.1) heeft. We kunnen de matrix-geometrische methode dus toepassen op het voorraadprobleem, waarbij we de toestanden opdelen in blokken ter grootte Q. We gebruiken de volgende notatie voor dit specifieke probleem: • H = {s + Q, s, s − Q, . . . }, • K = {0, 1, . . . , Q − 1}. Toestanden i in het voorraadmodel zijn dan van de vorm i = η − k, met η ∈ H en k ∈ K. Eveneens geven we de evenwichtsverdeling π ¯ als volgt weer: π ¯ = {¯ πs+Q , π ¯s , π ¯s−Q , . . . }, met π ¯η = (π(η, 0), π(η, 1), . . . , π(η, Q − 1)). Het bepalen van de evenwichtsverdeling π wordt dan gereduceerd tot het bepalen van π ¯s+Q , π ¯s en R, waarbij de andere π ¯η gegeven worden door π ¯s−iQ = π ¯ s Ri ,

i = 1, 2, . . . .

Eerst dienen we na te gaan wanneer de evenwichtsverdeling bestaat, oftewel wanneer het model positief recurrent is. Hiertoe gebruiken we de driftconditie (3.6). Herinner je

23

dat de matrices A0 , A1 en A2 gegeven worden door   0 0 ··· 0  .. .. ..   . A0 =  . . ,  0 0 · · · 0 λ 0 ··· 0  0 ··· 0 −λ − L1 λ ..  .. .. .. . . .  0 .  . . . .  .. .. .. A1 =  .. 0  . . . .. ..  .. λ 0 · · · · · · 0 −λ − 1  0 ··· 0 L . . ..  . ..  0 ..  A2 =  . . . . . . . . . . 0 0 · · · 0 L1



1 L

   ,   

Optellen van deze drie matrices leidt tot de uiterst eenvoudige matrix A = A0 + A1 + A2 :   −λ λ 0 ··· 0 ..   ... ... ...  0 .    . . . .  .. .. .. 0  A =  .. .   .. ..  0 . λ . λ 0 · · · 0 −λ Dit is simpelweg een cyclisch proces waarbij de tijd tussen transities overal exp(λ) verdeeld is. Het is gemakkelijk na te gaan dat de evenwichtsverdeling πA bij deze matrix gelijk is aan 1 1 πA = ,..., . Q Q De linkerkant van de driftconditie levert dus op πA A0 e =

λ . Q

De rechterkant van de ongelijkheid geeft Q X 1 1 πA A2 e = = . QL L i=1

Het proces is dus positief recurrent wanneer Qλ < L1 , oftewel wanneer λL < Q. Deze conditie is intu¨ıtief heel logisch, aangezien met deze conditie tekorten altijd afgelost kunnen worden door telkens te bestellen wanneer de vorige bestelling binnen is (dit gebeurt

24

wanneer de voorraad na binnenkomst van een bestelling onder s blijft). Immers, het verwachte aantal klanten in de levertijd is dan kleiner dan de bestelgrootte. We zullen er in het vervolg van uitgaan dat Q groot genoeg is zodat aan de driftconditie voldaan is. We richten ons nu op de berekening van de matrix R. We hebben gezien dat deze matrix aan de volgende kwadratische vergelijking moet voldoen: A0 + RA1 + R2 A2 = 0. Het direct bepalen van R uit deze vergelijking voor algemene Q is bijzonder lastig. Wij zullen R dus numeriek benaderen via de methode (3.10). We kunnen deze gebruiken aangezien A1 bovendriehoeks is met diagonaalelementen ongelijk aan nul (mits −λ − L1 6= 0), oftewel det(A1 ) 6= 0. Ter illustratie berekenen we numeriek de evenwichtsverdeling voor de volgende specifieke keuzes van de parameters: • λ = 2, • L = 1, • s = 3, • Q = 3. De matrixgeometrische methode groepeert de toestanden dan per drie, waarbij maximale voorraad 6 is. We hebben dus π ¯6 = (π6 , π5 , π4 ), π ¯3 = (π3 , π2 , π1 ), etc. De matrices A0 , A1 en A2 zijn dan als volgt:       0 0 0 −3 2 0 1 0 0 A0 = 0 0 0 , A1 =  0 −3 2  , A2 = 0 1 0 . 2 0 0 0 0 −3 0 0 1 Verder hebben we voor het bepalen van R de inverse van A1 nodig. Deze is in dit geval gelijk aan   1 4 − 3 − 29 − 27  0 −1 −2  . A−1 1 = 3 9 0 0 − 13 De matrices V en W in (3.10) zijn dan gelijk aan    1  4 0 0 0 − 3 − 29 − 27  0  0 −1 −2  . 0 0  , W = A2 A−1 V = A0 A−1 1 = 1 = 3 9 2 4 8 − 3 − 9 − 27 0 0 − 13

25

De eerste vier iteraties van R(k) zien er als volgt  0 0 0 R(1) =  0 0.6667 0.4444  0 0 0 R(2) =  0 0.7325 0.5322  0 0 0 R(3) =  0 0.7604 0.5751  0 0  0 0 R(4) = 0.7763 0.6005

uit:  0 0  0.2963  0 0  0.3841  0 0  0.4326  0 0 . 0.4627

We zien dat deze rij convergeert. Na veertig iteraties hebben we   0 0 0 0 0  R(40) =  0 0.8105 0.6570 0.5325 en bij verdere iteraties verschillen de elementen minder dan 0.001 van deze matrix. We zullen werken met deze matrix R. We moeten voor dit voorbeeld nog π ¯6 en π ¯3 berekenen. Via (3.11) zien we in dat deze de oplossing zijn van de volgende matrixvergelijking:   −2 2 0 0 0 0  0 −2 2  0 0 0     0 0 −2 2 0 0 B0 A0 . = (π6 , π5 , π4 , π3 , π2 , π1 )  (¯ π6 , π ¯3 ) 1  B1 A1 + RA2 0 0 −3 2 0   0  1 0 0 −3 2 0 0 1 0.8105 0.6570 −2.4675 Een oplossing van deze vergelijking is (ˆ π6 , π ˆ5 , π ˆ4 , π ˆ3 , π ˆ2 , π ˆ1 ) = (0.2311, 0.4184, 0.5703, 0.4622, 0.3747, 0.3037). Deze schalen we vervolgens met ˆ¯6 e + π ˆ¯3 (I − R)−1 e = 3.6596 ,α = π wat ons de oplossing van de eerste zes evenwichtskansen geeft: (π6 , π5 , π4 , π3 , π2 , π1 ) = (0.0632, 0.1143, 0.1558, 0.1263, 0.1024, 0.0830). Deze vector, samen met R, levert ons de volledige evenwichtsverdeling via de relatie π ¯3−3i = π ¯ 3 Ri ,

26

i = 1, 2, . . . .

De simpele vorm van R maakt het mogelijk de volgende expliciete uitdrukking voor de machten van de matrix te geven   0 0 0 0 0 0 , Ri =  i−1 i−1 0.8105 · 0.5325 0.6570 · 0.5325 0.5325i wat betekent dat we de andere πi ook zonder matrixvermenigvuldiging kunnen weergeven, en alleen gebruikmakend van π1 = 0.0830: π3−3i = 0.0830 · 0.8105 · 0.5325i−1 , i−1

π2−3i = 0.0830 · 0.6570 · 0.5325 π1−3i = 0.0830 · 0.5325i ,

,

i = 1, 2, . . . i = 1, 2, . . .

i = 1, 2, . . . .

Uiteindelijk kunnen we middels de evenwichtsverdeling eenvoudig de voorraadkosten en kosten bij tekorten berekenen. Deze kosten worden behandeld in de volgende sectie.

3.3. Voorraadkosten en kosten bij tekorten De evenwichtsverdeling van de voorraad levert ons vrijwel direct de voorraadkosten en kosten bij tekorten. We kunnen de evenwichtsverdeling namelijk opvatten als de fractie van tijd dat in elke toestand wordt besteed op lange termijn. De voorraadkosten en kosten bij tekorten zijn direct afhankelijk van deze fractie van tijd. Om precies te zijn, worden de gemiddelde voorraadkosten r¯ en gemiddelde boetekosten ¯b per tijdseenheid, gegeven voorraadkosten r per tijdseenheid per product in voorraad en boetekosten b per tijdseenheid per tekort, gedefinieerd door r¯ = lim E[fysieke voorraad op tijdstip t] · r

(3.12)

¯b = lim E[tekort op tijdstip t] · b

(3.13)

t→∞

t→∞

Maar voor deze verwachtingen hebben we enkel de evenwichtsverdeling nodig. Merk op dat de fysieke voorraad gelijk is aan de nettovoorraad wanneer de nettovoorraad positief is en 0 anders. Analoog geldt dat het tekort gelijk is aan -nettovoorraad als deze negatief is en 0 anders. Om de verwachtingen uit te rekenen gebruiken we zoals gebruikelijk een gewogen som, zodat we de volgende uitdrukkingen krijgen: lim E[fysieke voorraad op tijdstip t] =

t→∞

lim E[tekort op tijdstip t] =

t→∞

s+Q X i=1 ∞ X

iπi iπ−i

i=1

Beide sommen zijn in het geval van geometrisch gerelateerde evenwichtskansen gemakkelijk uit te rekenen. We zullen beide verwachtingen uitrekenen voor het voorbeeld dat we in de vorige sectie besproken hebben.

27

Eerst bepalen we de gemiddelde fysieke voorraad. De positieve evenwichtskansen waren de volgende: (π6 , π5 , π4 , π3 , π2 , π1 ) = (0.0632, 0.1143, 0.1558, 0.1263, 0.1024, 0.0830). P Wanneer we deze invullen in de eindige som, komen we uit op s+Q i=1 iπi = 2.2406. De gemiddelde voorraadkosten r¯ per tijdseenheid zijn in dit geval dus gelijk aan 2.2406r. Voor het gemiddelde tekort delen we de negatieve niveaus van de nettovoorraad op in drie stukken, ten einde de geometrische relatie te gebruiken die voor elke groep geldt. Er geldt: ∞ X

iπ−i =

i=1

∞ X

3iπ−3i +

i=1 ∞ X

=3

∞ X

(3i − 1)π1−3i +

i=1 ∞ X

iπ−3i + 3

i=1

∞ X

(3i − 2)π2−3i

i=1

iπ1−3i + 3

i=1

∞ X

iπ2−3i −

i=1

∞ X

π1−3i − 2

i=1

∞ X

π2−3i

i=1

Al deze oneindige sommen hebben een van de volgende vormen: ∞ X i=1

a , ia = (1 − a)2 i

∞ X i=1

ai =

a . 1−a

Deze gelijkheden kan men eenvoudigweg nagaan door de sommen met ofwel (1 − a)2 of 1 − a te vermenigvuldigen. Door deze identiteiten toe te passen op deze vijf sommen, verkrijgen we uiteindelijk de volgende waarde: ∞ X

iπ−i = 0.4917 + 0.6067 + 0.7485 − 0.0945 − 0.2333 = 1.5191.

i=1

De gemiddelde boetekosten zijn dus b = 1.5191b in dit voorbeeld. Deze berekeningen kunnen worden uitgevoerd voor een willekeurige combinatie van s, Q, λ en L. De voorraadkosten en kosten bij tekorten zijn middels de evenwichtskansen dus relatief gemakkelijk te berekenen. Het bepalen van het geld dat gemiddeld wordt besteed aan de vaste kosten per bestelling vereist een iets andere aanpak.

3.4. Vaste kosten De berekening van de vaste kosten op de lange termijn verloopt anders dan die van de andere kosten, doordat vaste kosten niet verbonden zijn aan een toestand, maar aan een toestandsovergang. Specifiek worden vaste kosten opgelopen wanneer het proces van toestand i naar i + Q gaat. De bijbehorende overgangsintensiteit wordt gegeven door qi,i+Q . Deze overgangsintensiteit geeft het gemiddelde aantal overgangen naar toestand i + Q aan, startend vanuit toestand i. Hierdoor worden de gemiddelde vaste kosten per toestand i gegeven door qi,i+Q K. Op de lange termijn besteedt het proces een fractie πi van de tijd in toestand i. De vaste kosten worden nu bepaald door de gemiddelde kosten in elke toestand te schalen

28

met de evenwichtsverdeling. De gemiddelde kosten per tijdseenheid worden dus gegeven door: ! Q X X X X K K 1 ¯ = πs+i . (3.14) πi = 1− K πi qij cij = πi · · K = L L i≤s L i=1 i≤s i≤s+Q Deze kosten zijn computationeel gemakkelijk uit te rekenen. We hoeven enkel een eindige som uit te rekenen met evenwichtskansen die we hebben bepaald middels de matrixgeometrische methode. Merk verder op dat de gemiddelde vaste kosten onafhankelijk zijn van het bestelpunt s, aangezien de evenwichtskansen πs+i voor i ∈ {1, . . . , Q} onafhankelijk zijn van s. We berekenen ook de vaste kosten voor het voorbeeld s = Q = 3. De eindige som is in dit geval gelijk aan π4 + π5 + π6 = 0.3333. De gemiddelde vaste kosten zijn dus ¯ = K (1 − 0.3333) = 0.6667K. K L

3.5. Combinaties van s en Q We hebben nu methoden gevonden om numeriek de verscheidene types kosten te berekenen voor willekeurige s, Q, λ en L. Het optimaliseren van deze kosten voor het (s, Q)-voorraadmodel is, gezien het feit dat deze numeriek bepaald zijn en niet analytisch, buiten het bereik van dit onderzoek. Wel zal het inzichtelijk zijn deze kosten weer te geven voor verschillende combinaties van s en Q. Hierbij stellen we de volgende parameters vast: • λ = 2, • L = 1, • r = 1, • K = 1, • b = 3. Voor deze parameters bepalen we de totale gemiddelde kosten per tijdseenheid, en wel voor alle combinaties van s ∈ {1, . . . , 5} en Q ∈ {3, . . . , 7} (Q > 2 is nodig voor positieve recurrentie). Hiertoe voeren we hetzelfde procedé uit dat we voor het voorbeeld met s = Q = 3 gebruikt hebben. We bepalen eerst voor elke combinatie de gemiddelde voorraadkosten:

29

1 2 3 4 5

s

3 1.0336 1.5956 2.2406 2.9529 3.7197

4 1.6969 2.4222 3.2186 4.0677 4.9559

Q 5 2.2548 3.0538 3.9113 4.8103 5.7387

6 2.7712 3.6115 4.5010 5.4245 6.3715

7 3.2706 4.1375 5.0467 5.9847 6.9424

Tabel 3.1.: Voorraadkosten r¯ We zien dat aanpassingen van s grotere effecten hebben op de voorraadkosten, naarmate s groter wordt. Daarentegen hebben aanpassingen van Q steeds kleinere effecten op de voorraadkosten bij grotere Q. Verder zorgt een hoger bestelpunt voor grotere invloeden van aanpassingen van Q op de prijs. Dit fenomeen is vooral aanwezig bij lage bestelgroottes. Verder bekijken we de gemiddelde kosten bij tekorten, in de vorm van boetekosten. Weer geven we deze voor verschillende combinaties van s en Q weer.

1 2 3 4 5

s

3 6.9367 5.6225 4.5573 3.7410 3.2901

4 3.1859 2.3616 1.7505 1.2976 1.0027

Q 5 2.0725 1.4693 1.0417 0.7386 0.5237

6 1.5549 1.0760 0.7447 0.5155 0.3567

7 1.2589 0.8594 0.5868 0.4006 0.2735

Tabel 3.2.: Boetekosten ¯b We zien dat voor de kleinere bestelgroottes, het vergroten van de bestelgrootte de boetekosten behoorlijk inperkt. Dit effect is ook terug te zien in het bestelpunt, zij het in mindere mate. Verder heeft het verhogen van het bestelpunt als gevolg dat aanpassingen van de bestelgrootte minder effect heeft op de boetekosten. Evenzo hebben verhogingen van het bestelpunt minder effect wanneer de bestelgrootte toeneemt. ¯ Aangezien deze Ten slotte bekijken we de gemiddelde vaste kosten per tijdseenheid K. onafhankelijk zijn van het bestelpunt, hoeven we slechts vijf waarden te berekenen. De gemiddelde vaste kosten voor elke bestelhoeveelheid worden weergegeven in de volgende tabel: Q ¯ K

3 4 5 6 7 0.6667 0.5000 0.4000 0.3333 0.2857 ¯ Tabel 3.3.: Vaste kosten K

We zien dat het verhogen van de bestelgrootte een verminderd effect heeft op de vaste kosten naarmate Q groter wordt. In dit voorbeeld specifiek zien we een patroon in de

30

¯ = 2 . De auteur verwacht dat de vaste kosten. Ze worden in dit geval gegeven door K Q gemiddelde vaste kosten per tijdseenheid gelijk zullen zijn aan ¯ = Kλ . K LQ Het bewijzen van deze formule is een onderwerp voor verder onderzoek. Uiteraard zijn we uiteindelijk ge¨ınteresseerd in de gemiddelde totale kosten per tijdseenheid. Om in de totale kosten inzicht te krijgen combineren we de drie tabellen. De volgende tabel geeft de totale gemiddelde kosten aan voor de combinaties van s en Q:

s

1 2 3 4 5

3 7.6370 7.8848 7.4646 7.3610 7.6765

4 5.3828 5.2838 5.4691 5.8653 6.4586

Q 5 4.7273 4.9231 5.3530 5.9489 6.6624

6 4.7261 5.0208 5.5790 6.2733 7.0615

7 4.8152 5.2826 5.9192 6.6710 7.5016

Tabel 3.4.: Totale kosten Het blijkt zeer onvoordelig te zijn een waarde van Q te kiezen die net boven λL ligt (in dit geval Q = 3). Dit komt doordat de voorraad dan onnodig lang leeg is en het lang duurt voordat tekorten gecompenseerd worden. Relatief hoge waarden van Q zijn wenselijk. Verder zien we dat de kosten sterker in s variëren wanneer Q groter is. Het optimum voor s en Q lijkt te liggen bij s = 1 en Q = 6. Er is geen aanleiding te vinden dat het optimum ergens anders ligt. Immers, s = 1 is het minimale bestelpunt dat gekozen kan worden in dit model, en hogere bestelpunten dan s = 1 leiden in dit geval tot hogere kosten, behalve in het geval Q = 4. Ook voor de bestelgrootte geldt dat bestelgroottes die niet gelijk zijn aan 6 leiden tot hogere kosten, tenzij het bestelpunt hoger wordt. We hebben nu inzicht gekregen in de structuur van de gemiddelde kosten in ons model. We hebben zelfs een optimale (s, Q)-strategie gevonden voor specifieke parameterkeuzes. Maar dit optimum is slechts gevonden door het uitproberen van waarden voor s en Q. Bovendien: de optimale bestelstrategie hoeft helemaal niet een (s, Q)-strategie te zijn. Immers, bij de introductie van ons model hebben we al aangekaart dat waarschijnlijk wenselijk is de bestelgrootte niet vast te zetten. Aangezien er één bestelling tegelijk kan worden geplaatst, kan het zijn dat de voorraad onder s gedaald is op het moment dat de volgende bestelling weer mogelijk is. Om tekorten te beperken, zal er dan meer besteld moeten worden. In het volgende hoofdstuk laten we de restrictie los dat de bestelgrootte vaststaat. Dit stelt ons in staat ons voorraadmodel te formuleren als een zogenaamd Markov beslissingsproces. Het vakgebied dat deze processen bestudeert, de Markov beslissingstheorie, heeft als doel een optimale strategie te vinden voor zulke processen. Hiertoe zijn verscheidene algoritmes ontwikkeld, die numeriek geimplementeerd kunnen worden. We behandelen een dergelijk algoritme dat ons in staat stelt het voorraadprobleem op te lossen.

31

4. Markov beslissingsproces voor voorraad Markov Decision Processes, of Markov beslissingsprocessen zijn zeer populaire manieren om een stochastisch optimalisatieprobleem te modelleren. In het geval van stochastische processen gedreven door Markovketens, zoals het probleem dat wij bekijken, blijkt deze theorie een veelzijdige aanpak te geven. We zullen het voorraadproces modelleren via een MDP om vervolgens de technieken toe te passen die ontwikkeld zijn voor dit soort processen. Voordat we een MDP kunnen construeren, bouwen we de theorie op.

4.1. Markov beslissingstheorie Een MDP wordt als volgt gedefinieerd. Definitie 4.1. Een MDP bestaat uit een (discrete-tijds) proces Xn , n = 0, 1, 2, . . . met toestandsruimte I en een reeks acties An , n = 0, 1, 2, . . . met An ∈ A die het verloop van Xn be¨ınvloedt. De verzameling A heet de actieruimte. Aanname 4.2. De toestandsruimte I en de actieruimte A zijn aftelbaar. Vaak zal het geval zijn dat er in een bepaalde toestand Xn = i slechts een deelverzameling van acties An ∈ Ai ⊂ A toegestaan zijn. Verder nemen we de volgende welbekende eigenschap aan. Definitie 4.3. De processen Xn , An , n = 0, 1, 2, . . . in een MDP voldoen aan de volgende eigenschap: P(Xn+1 = in+1 |X0 = i0 , . . . , Xn = in , A0 = a0 , . . . , An = an ) = P(Xn+1 = in+1 |Xn = in , An = an ).

(4.1) (4.2)

Net als in een normale Markovketen hangt de volgende toestand Xn+1 dus alleen af van de huidige toestand van Xn en An en niet van het verleden. Deze aanname roept de volgende notatie op: pa (i, j) := P(Xn+1 = j|Xn = i, An = a).

(4.3)

De n-stapskansen worden weergegeven met pan (i, j). Uiteraard geldt pa1 (i, j) = pa (i, j). Aan elke actie a zijn zekere beloningen ra (i, j) verbonden. Definieer daarmee:

32

Definitie 4.4. Stel dat een beloning ra (i, j) verdiend wordt wanneer het proces Xn zich in toestand i bevindt, actie a genomen wordt en het proces naar toestand j gaat. Dan definieert X ra (i) := pa (i, j)ra (i, j) (4.4) j∈I

de verwachte beloning als actie a wordt genomen in toestand i. Aanname 4.5. De verwachte beloningen zijn uniform begrensd, dus er bestaat een R > 0 zodanig dat ||ra (i)|| < R voor alle i ∈ I en a ∈ Ai . Het probleem dat vaak opgelost moet worden in een MDP is de maximalisatie van de verwachte beloningen na een bepaalde tijd T (of de gemiddelde beloning per tijdseenheid voor T → ∞). Verder is er vaak sprake van een eindbeloning q(i) voor wanneer XT = i, die meegenomen wordt in het probleem. Voor deze maximalisatie worden strategieën gebruikt. Een strategie definieert een actie die genomen wordt voor elke toestand waarin het systeem zich bevindt. Ze worden gegeven door functies die aan elke toestand i een actie a toekennen. De volgende definitie formaliseert dit begrip. Definitie 4.6. Een strategie is een functie f die de toestandsruimte I afbeeldt op de actieruimte A. De gekozen strategie heeft invloed op het verloop van het stochastische proces Xn en op de beloningen. We zullen de volgende verkorte notatie gebruiken voor de overgangskansen en beloningen: pf (i, j) := pf (i) (i, j), rf (i) = rf (i) (i), i, j ∈ I. De uiteindelijke grootheid die gemaximaliseerd dient te worden is de verwachte beloning tot tijdstip T gegeven een strategie f en startend vanaf tijdstip 0: VTf (i)

:=

T −1 X

Ef [rAn (Xn , Xn+1 )|X0 = i] + Ef [q(XT )|X0 = i].

(4.5)

n=0

Merk op dat de verwachting van de beloningen afhankelijk is van de gebruikte strategie. Een superscript f geeft deze afhankelijkheid aan. Middels de definitie van de verwachte beloning (4.5) kunnen we VTf ook als volgt schrijven: VTf (i)

=

T −1 X

Ef [rAn (Xn )|X0 = i] + Ef [q(XT )|X0 = i].

(4.6)

n=0

Als we het probleem bekijken met oneindige tijdshorizon, dus als T → ∞, dan maximaliseren we de gemiddelde verwachte beloning g f (i): g f (i) := lim sup T →∞

T −1 1 X f An VTf (i) = lim sup E [r (Xn )|X0 = i], T T →∞ T n=0

waarbij we stellen dat de eindbeloningen q(i) gelijk zijn aan 0 voor alle i ∈ I.

33

(4.7)

We noteren VT∗ (i) en g ∗ (i) als de optimale waarden van VTf (i) en g f (i) respectievelijk. Dat wil zeggen, (4.8) VT∗ (i) := sup VTf (i), g ∗ (i) := sup g f (i). f

f

Wij zullen voor het optimaliseren van onze voorraadstrategie kijken naar een oneindige tijdshorizon, aangezien we de kosten willen minimaliseren op de lange termijn. Een nuttige methode om de optimale strategie te benaderen maakt niet gebruik van de evenwichtskansen maar maakt gebruik van een iteratief algoritme. De Markov-beslissingstheorie levert een aantal krachtige algoritmes om een optimale strategie te bepalen. Wij zullen gebruik maken van de volgende methode, die in de volgende sectie theoretisch besproken wordt.

4.2. Successieve approximatie De successieve-approximatiemethode is gebaseerd op een zeer krachtige stelling in de theorie van Markov-beslissingsprocessen. Deze stelling geeft een vergelijking die een optimale strategie kenmerkt. Stelling 4.7. Stel er bestaat een begrensde functie d(i), i ∈ I, en een constante g zodanig dat voor alle i ∈ I geldt, ( ) X d(i) + g = max ra (i) + pa (i, j)d(j) . (4.9) a∈Ai

j∈I

Dan g = g ∗ (i) voor alle i ∈ I en dan geldt dat elke strategie f , die voldoet aan, ( ) X f (i) ∈ arg max ra (i) + pa (i, j)d(j) , a∈Ai

j∈I

de gemiddelde beloning maximaliseert: g f (i) = g ∗ . Bewijs. Het bewijs van deze stelling wordt toegeschreven aan R. N´ un ˜ez-Queija[2]. We V f (i)

bewijzen eerst dat de gemiddelde beloning na lim supT →∞ TT begrensd wordt door g en daarna dat deze begrenzing bereikt kan worden. We bekijken d(·) als functie van Xt en bepalen het verwachte verschil van d(Xt+1 ) met d(Xt ) gegeven de waarde van Xt . Voor elke strategie f geldt X Ef [d(Xt+1 )|Xt = i] = pf (i, j)d(j) + rf (i) − rf (i) j

! ≤ max a

X

a

p (i, j)d(j) + r (i)

j

= d(i) + g − rf (i),

34

a

− rf (i)

voor alle i ∈ I. De verwachte stijging van d(·) na één stap heeft dus als bovengrens g − rf (i). Middels deze bovengrens kunnen we een bovengrens vinden voor de stijging van d(·) in de (t + 1)-de stap, door te conditioneren op de toestand op tijdstip t: X f Ef [d(Xt+1 )|X0 = i] = pt (i, j)Ef [d(Xt+1 )|Xt = j] j

≤

X

pft (i, j)(d(j) + g − rf (i))

j

= g + Ef [d(Xt )|X0 = i] − Ef [rAt (Xt )|X0 = i], voor alle i ∈ I. Als we nu deze ongelijkheid sommeren van t = 0 tot T − 1 verkrijgen we T −1 X

f

E [d(Xt+1 )|X0 = i] ≤ T g +

t=0

= Tg +

T −1 X t=0 T −1 X

f

E [d(Xt )|X0 = i] +

T −1 X

Ef [rAt (Xt )|X0 = i]

t=0

Ef [d(Xt )|X0 = i] + VTf (i),

t=0

waaruit volgt Ef [d(XT )|X0 = i] ≤ T g + Ef [d(X0 )|X0 = i] + VTf (i), en aangezien d(i), i ∈ I een begrensde functie is, geldt dus lim sup T →∞

VTf (i) ≤ g. T

Het bewijs van de existentie van een strategie f die deze bovengrens bereikt, vereist enkel de opmerking dat alle ongelijkheiden in de bovenstaande vergelijkingen gelijkheden worden wanneer de strategie f voldoet aan ( ) X f (i) ∈ arg max ra (i) + pa (i, j)d(j) . a∈Ai

j∈I

Het successieve-approximatie (SA) algoritme maakt gebruik van deze conditie door de functie d(i) op te bouwen in een groot aantal stappen en telkens de vorige iteratie te gebruiken in de huidige. Het algoritme werkt als volgt. Het algoritme 0. Laat n := 0. Kies een > 0 en een begrensde functie v0 (i) (vaak wordt v0 (i) ≡ 0 gebruikt). 1. Bereken

(

) a

vn+1 (i) := max r (i) + a∈Ai

X j∈I

35

a

p (i, j)vn (j)

(4.10)

en laat

( fn+1 (i) ∈ arg max ra (i) + a∈Ai

) X

pa (i, j)vn (j) .

(4.11)

j∈I

2. Laat Mn := maxi∈I {vn (i) − vn−1 (i)} en mn := mini∈I {vn (i) − vn−1 (i)}. Stop het algoritme als Mn − mn < . Anders laten we n := n + 1 en herhalen we stappen 1 en 2. Als we v0 (i) ≡ 0 kiezen, is het intu¨ıtief niet lastig in te zien dat dit algoritme de optimale strategie benadert. Immers, vn (i) kan dan worden ge¨ınterpreteerd als de maximale beloning over n perioden, waardoor de maximale gemiddelde beloning per tijdsstap benaderd wordt door vn (i) − vn−1 (i). Om te laten zien dat dit algoritme daadwerkelijk de optimale strategie benadert, doen we eerst een aanname met betrekking tot de strategieën. Deze aanname heeft te maken met de verwachte aankomsttijden T f (i, i0 ) bij een vaste toestand i0 vanaf toestand i, afhankelijk van de gebruikte strategie f : T f (i, i0 ) = Ef [inf{n ≥ 1 : Xn = i0 }|X0 = i]. Aanname 4.8. De Markovketen met overgangskansen pf (i, j) is aperiodiek. Verder kan er een vaste toestand i0 en T0f < ∞ gekozen worden zodanig dat T f (i, i0 ) < T0f voor alle i ∈ I. Hieruit volgt dat de toestand i0 positief recurrent is. Het volgende lemma helpt ons te bewijzen dat het SA-algoritme convergeert naar de optimale strategie. Lemma 4.9. Laat f een strategie zijn die voldoet aan Aanname (4.8) en stel dat voor een zekere constante g en begrensde functie v(i), i ∈ I geldt dat X rf (i) + pf (i, j)v(j) ≥ v(i) + g, i ∈ I. (4.12) j∈I

Dan g f ≥ g. Analoog, als rf (i) +

X

pf (i, j)v(j) ≤ v(i) + g,

i ∈ I,

(4.13)

j∈I

dan g f ≤ g. Bewijs. Het bewijs van dit lemma, toegewezen aan H.C. Tijms[3], gaat analoog voor beide richtingen van de ongelijkheid. Bekijk dus het geval dat er voldaan is aan (4.12). We bewijzen de volgende bewering middels inductie naar T : X f VTf (i) − T g + pT (i, j)v(j) ≥ v(i), i ∈ I, T = 1, 2, . . . . (4.14) j∈I

36

De basisstap is triviaal. Voor T = 1 komen de ongelijkheden (4.12) en (4.14) immers precies overeen. Voor T ≥ 2 conditioneren we de overgangskansen pfT (i, j) op pfT −1 (i, j). Het is gemakkelijk in te zien dat X f pfT (i, j) = pT −1 (i, k)pf (k, j). k∈I

Verder gebruiken we de definitie van de verwachte beloningen tot op tijdstip T: VTf (i)

=

T −1 X

Ef [rAn (Xn )|X0 = i]

n=0

Volgens de inductiehypothese geldt de bewering voorP T − 1. Als we nu rf (i) in (4.12) naar rechts halen en de ongelijkheid substitueren voor j pf (i, j)v(j), verkrijgen we VTf (i) − T g +

X

pfT (i, j)v(j)

j∈I

= VTf (i) − T g +

XX

pfT −1 (i, k)pf (k, j)v(j)

j∈I k∈I

=

VTf (i)

− Tg +

X

pfT −1 (i, k)

≥ VTf (i) − T g +

X

pf (k, j)v(j)

j∈I

k∈I 4.12

X

pfT −1 (i, k)(v(k) + g − rf (k))

k∈I

=

VTf (i)

−

X

pfT −1 (i, k)rf (k)

− (T − 1)g +

k∈I

X

pfT −1 (i, k)v(k)

k∈I

= VTf (i) − Ef [rAT −1 (XT −1 )|X0 = i] − (T − 1)g +

X

pfT −1 (i, k)v(k)

k∈I

= VTf−1 (i) − (T − 1)g +

X

IH

pfT −1 (i, k)v(k) ≥ v(i).

k∈I

Hiermee is 4.14 bewezen. Het bewijs van het lemma volgt nu vrijwel direct: deel beide kanten van (4.14) door T en laat T → ∞. Er volgt dan dat g f − g ≥ 0, wat de eerste helft van het lemma bewijst. De tweede helft volgt analoog door de ongelijkheden in het bewijs om te draaien. De volgende stelling geeft begrenzingen voor de gemiddelde beloning die de strategieën fn opleveren die in het SA-algoritme gegenereerd worden. Stelling 4.10. Laat vn (i), fn (i), Mn en mn berekend zijn uit het SA-algoritme. Als fn voldoet aan Aanname (4.8) dan mn ≤ g fn ≤ g ∗ ≤ Mn , en mn is niet-dalend en Mn niet-stijgend.

37

Bewijs. Ook het bewijs van deze stelling wordt geattribueerd aan Tijms[3]. Vanwege Lemma (4.9) is het voldoende de ongelijkheid X rfn (i) + pfn (i, j)vn−1 (j) ≥ vn−1 (i) + mn , i ∈ I, (4.15) j∈I

te bewijzen om te laten zien dat mn ≤ g fn . Voor Mn bewijzen we X rf (i) + pf (i, j)vn−1 (j) ≤ vn−1 (i) + Mn , i ∈ I,

(4.16)

j∈I

voor alle strategieën f . Immers, dan geldt dat g f ≤ Mn voor alle strategieën f en dus dat g ∗ ≤ Mn . Om de eerste ongelijkheid te bewijzen, merken we op dat uit de definitie van fn (i) volgt dat X vn (i) = rfn (i) + pfn (i, j)vn−1 (j), i ∈ I. (4.17) j∈I

Door nu te gebruiken dat vn (i) − vn−1 (i) ≥ mn voor alle i verkrijgen we X rfn (i) + pfn (i, j)vn−1 (j) = vn (i) = vn−1 (i) + vn (i) − vn−1 (i) ≥ vn−1 (i) + mn , j∈I

wat precies de ongelijkheid is die we zochten. Voor de tweede ongelijkheid merken we op dat vn in het SA-algoritme zo gekozen is dat voor alle strategieën f geldt, X pf (i, j)vn−1 (j) ≤ vn (i), i ∈ I. (4.18) rf (i) + j∈I

Weer vullen we vn (i) = vn−1 (i)+vn (i)−vn−1 (i) in en nu gebruiken we dat vn (i)−vn−1 (i) ≤ Mn voor alle i om de gewenste ongelijkheid te verkrijgen. We moeten nu alleen nog bewijzen dat mk+1 ≥ mk en Mk+1 ≤ Mk . Laat n = k in (4.17) en n = k + 1 en f = fk in (4.18) en trek de vergelijkingen van elkaar af. Dan geldt X vk+1 (i) − vk (i) ≥ pfk (i, j)(vk (j) − vk−1 (j)), i ∈ I. j∈I

Hieruit volgt dat vk+1 (i) − vk (i) ≥ mk voor alle i ∈ I en dus dat mk+1 ≥ mk . Andersom, laat n = k + 1 in (4.17) en n = k en f = fk+1 in (4.18). Het verschil van de vergelijkingen levert dan op X vk+1 (i) − vk (i) ≤ pfk+1 (i, j)(vk (j) − vk−1 (j)). j∈I

Evenzo volgt hieruit dat vk+1 (i) − vk (i) ≤ Mk voor alle i ∈ I en dus dat Mk+1 ≤ Mk . Dit betekent dat de gemiddelde beloning g fn convergeert naar het optimum g ∗ als Mn − mn convergeert naar 0 en dat, in dat geval, het optimum tevens benaderd wordt door mn en Mn . We zullen zien dat dit in ons voorraadmodel het geval is.

38

Om dit algoritme te gebruiken heeft de auteur code geschreven in Matlab dat de verscheidene berekeningen uitvoert en stopt wanneer het gewenste nauwkeurigheidsniveau is bereikt (wanneer Mn − mn < ). We gaan er hierbij eerst van uit dat de actieruimte eindig is. Dit is nodig zodat een computer de maxima in (4.10) en (4.11) kan berekenen. Nu we een techniek hebben om een optimale strategie te bepalen, zullen we het voorraadprobleem formuleren als een MDP om vervolgens het SA-algoritme toe te passen.

4.3. Formulatie MDP De eerste taak die we moeten vervullen voor de formulatie van het voorraadprobleem als MDP is een beschrijving te geven voor het stochastische voorraadproces als Markovketen. Aangezien de levertijden en tijden tussen aankomsten van klanten beiden exponentieel verdeeld zijn in dit model, kunnen we het proces beschrijven als een continue-tijds Markovketen. We kunnen er echter niet meer in volstaan alleen de nettovoorraad bij te houden, aangezien we de bestelgrootte Q niet meer vast houden. We voegen een tweede dimensie toe aan de toestand van het proces, namelijk de bestelgrootte van de te leveren bestelling (deze is 0 wanneer er geen bestelling in behandeling is). Een toestand heeft dan de vorm (i, x), waarbij i ∈ Z de nettovoorraad is, en x ≥ 0 de bestelgrootte van de huidige bestelling. Uiteraard moeten ook de genomen acties a ∈ A := Z≥0 bijgehouden worden, aangezien deze het voorraadproces be¨ınvloeden. De overgangsintensiteiten hebben dus vijf argua menten en zijn van de vorm q(i,x),(j,y) , waarbij i en j voorraadsniveaus aangeven, x en y bestelde hoeveelheden, en a de huidige actie (hoeveel bestellen). Dit lijkt veel om bij te houden, maar de mogelijke acties hangen zeer nauw samen met de bestelgrootte x van de huidige bestelling. Zo kan er geen a > 0 gekozen worden wanneer x > 0, aangezien er maar één bestelling tegelijk geplaatst kan worden. Het gedrag van het resulterende proces is als volgt. • Gedurende het gehele voorraadproces loopt een sterfteproces dat de nettovoorraad i met 1 doet afnemen met intensiteit λ (aankomst van klanten). • Wanneer de bestelde hoeveelheid x positief is, loopt er verder een proces dat x optelt bij de nettovoorraad na een levertijd die exp( L1 ) verdeeld is. Dit proces zorgt tevens ervoor dat de bestelde hoeveelheid 0 wordt. • Er kunnen acties (bestellingen) gedaan worden bij elke sprong die het proces maakt. Voor x > 0 is alleen de actie a = 0 mogelijk (niets bestellen). Als x = 0, kan er een hoeveelheid a > 0 besteld worden. • Voor acties a > 0 loopt eveneens een proces dat x optelt bij de nettovoorraad na een exp( L1 ) verdeelde levertijd. Als een actie a > 0 gekozen is, wordt er x := a gesteld bij de volgende sprong, mits een klant arriveert voor aankomst van de bestelling. Anders blijft x = 0 en kan er weer een actie a > 0 gedaan worden.

39

In termen van de overgangsintensiteiten hebben we dus 0 q(i,x),(i−1,x) = λ, 1 0 q(i,x),(i+x,0) = , L a q(i,0),(i−1,a) = λ, 1 a q(i,0),(i+a,0) = , L 0 q(i,0),(i,0) = −λ, a 0 q(i,0),(i,0) = q(i,x),(i,x)

i ∈ Z, x ≥ 0, i ∈ Z, x > 0, i ∈ Z, a > 0, i ∈ Z, a > 0,

i ∈ Z, 1 = −λ − , i ∈ Z, x > 0, a > 0. L

a waarbij de andere q(i,x),(j,y) met i 6= j en x 6= y gelijk zijn aan 0. Om het proces dan compleet te definiëren, geven we de mogelijke acties per toestand aan met de verzamelingen A(i,x) . Aangezien er één bestelling tegelijk geplaatst kan worden, worden deze simpelweg gegeven door

A(i,0) = A = Z≥0 , i ∈ Z, A(i,x) = {0}, i ∈ Z, x > 0. In de theorie voor MDP’s die we in de vorige sectie besproken hebben, zijn we uitgegaan van stochastische processen in discrete tijd. Het voorraadprobleem hoort bij een continuetijds Markovketen, maar toch kunnen we deze beschouwen als een discrete keten, door middel van uniformisatie.

4.3.1. Uniformisatie Laat Q = (qij )i,j∈I de genererende matrix zijn van een continue-tijds Markovketen.. Als de overgangsintensiteiten qii uniform begrensd zijn door een zekere constante γ, dat wil zeggen, γ ≥ sup |qii |, i

dan heeft een CTMC tevens een andere interpretatie. We kunnen een Markovketen dan opvatten als een keten waarbij we ’sprongen’ na dezelfde toestand toelaten, en de tijd tussen sprongen overal exp(γ) verdeeld is. De overgangskansen (pij ) worden dan gegeven door (q ij als i 6= j, pij = γ P qik 1 − k6=i γ als i = j. Dit proces heet uniformisatie. Aangezien de overgangsintensiteiten in het voorraadprobleem uniform begrensd zijn door γ := λ + L1 , kunnen we derhalve het voorraadproces

40

beschrijven door de volgende overgangskansen (pa(i,x),(j,y) )i,j∈Z,x,y,a≥0 : p0(i,x),(i−1,x) = p0(i,x),(i+x,0) pa(i,0),(i−1,a) pa(i,0),(i+a,0)

=

λ λ+ 1 L

λ+ λ = λ+ =

1 L 1 L 1 L

1 L

λ+

1 L

=

λL , λL + 1

1 , λL + 1 λL , = λL + 1 =

=

1 , λL + 1

i ∈ Z, x ≥ 0, i ∈ Z, x > 0, i ∈ Z, a > 0, i ∈ Z, a > 0,

λL 1 = , i ∈ Z, λL + 1 λL + 1 λL 1 =1− − = 0, i ∈ Z, x > 0, a > 0, λL + 1 λL + 1 = 0, elders.

p0(i,0),(i,0) = 1 − pa(i,0),(i,0) = p0(i,x),(i,x) pa(i,x),(j,y)

Hiermee hebben we de overgangskansen in de MDP-formulatie van ons probleem volledig beschreven. We moeten nu enkel nog de kosten ra (i, x) definiëren die verbonden zijn aan elke toestand en bijbehorende actie, voordat we het probleem aanpakken met het SA-algoritme.

4.3.2. Kosten in MDP We hoeven voor de formulatie van de kosten in een MDP enkel de verwachte kosten in elke toestand te berekenen. Aangezien het voorraadmodel een continue-tijds Markovketen is (die we discreet kunnen maken middels uniformisatie), zijn deze kosten enkel afhankelijk van de huidige toestand van het model. We moeten uitdrukkingen vinden voor de volgende drie soorten kosten: • Voorraadkosten rra (i, x) (r per tijdseenheid per eenheid product in voorraad) a • Vaste kosten per bestelling rK (i, x) (K per bestelling)

• Boetekosten rba (i, x) (b per tijdseenheid per eenheid product tekort) Eerst bekijken we de voorraadkosten en boetekosten, aangezien deze op een analoge manier berekend worden. De tijd tussen sprongen in de ge¨ uniformiseerde keten is voor elke L toestand exp(λ + L1 ) verdeeld. De verwachte tijd tussen overgangen is dus λ+1 1 = λL+1 . L De voorraadkosten rra (i, x) en boetekosten rba (i, x) worden dan simpelweg gegeven door ( L · ir als i > 0 a λL+1 , x ≥ 0, a ≥ 0, rr (i, x) = 0 als i ≤ 0 ( L − λL+1 · ib als i < 0 a rb (i, x) = , x ≥ 0, a ≥ 0. 0 als i ≥ 0

41

Verder hebben we nog vaste kosten per bestelling. In het MDP hebben we een gemakkelijke manier om bestellingen te identificeren. Er worden namelijk vaste kosten betaald telkens wanneer a > 0. Er geldt dus a rK (i, x) = K,

i ∈ Z, x ≥ 0, a > 0.

Combineren we de verschillende kosten, dan hebben we de volgende uitdrukking voor ra (i, x):  L · ir als i ≥ 0, a = 0,  λL+1    L · ir + K als i ≥ 0, a > 0, , x ≥ 0. (4.19) ra (i, x) = λL+1L  · ib als i < 0, a = 0, −  λL+1   L · ib + K als i < 0, a > 0, − λL+1 Merk op dat we hier te maken hebben met kosten. De Markov beslissingstheorie ging uit van beloningen die gemaximaliseerd dienden te worden. In dit geval moeten de kosten geminimaliseerd worden. Door de kosten te interpreteren als negatieve beloningen zien we echter in dat we in dit geval simpelweg de maxima in het SA-algoritme kunnen vervangen door minima. Het algoritme blijft verder onveranderd en benadert nog steeds de optimale strategie. Hiermee hebben we alle parameters gedefinieerd die we nodig hebben om het voorraadmodel als MDP te behandelen. Om daadwerkelijk het SA-algoritme toe te passen op het voorraadprobleem met behulp van de computer, is het echter nog nodig een maximum in te stellen op de bestelgrootte en de toestandsruimte te begrenzen. De stappen die we hiertoe moeten nemen worden besproken in de volgende sectie.

4.4. Truncatie van de toestandsruimte Om het SA-algoritme toe te passen middels de computer, hebben we een eindige actieruimte en toestandsruimte nodig. Om de toestandsruimte te begrenzen, moeten we echter ook de overgangskansen en mogelijke acties per toestand aanpassen. Eerst bekijken we de overgangskansen pa(i,x),(j,y) . Stel dat we een bovengrens M en een ondergrens m instellen voor de mogelijke niveaus van de nettovoorraad. Voor i = m kunnen we dan niet meer een positieve kans toelaten dat er een klant arriveert die de voorraad doet dalen. Om toch de kansen tot 1 te laten opsommen, moeten we de kansen om op voorraad m te blijven aanpassen. We hebben dan: p0(m,0),(m,0) = 1, λL , a > 0, x > 0. λL + 1 Om de bovengrens te realiseren, moeten we de mogelijke bestelgroottes beperken. Aangezien we niet boven een voorraad van M mogen komen, mag er voor een huidige voorraad van i niet meer dan M − i worden besteld. In termen van de actieruimtes A(i,x) hebben we dus pa(m,0),(m,0) = p0(m,x),(m,x) =

A(i,0) = {0, 1, 2, . . . , M − i}, m ≤ i ≤ M, A(i,x) = {0}, m ≤ i ≤ M, x > 0.

42

Uiteraard geldt hierdoor voor de bestelde hoeveelheden x een soortgelijke restrictie, waardoor ook de bestelde hoeveelheden een eindige verzameling vormen. De toestandsruimte en actieruimte zijn nu effectief eindige verzamelingen van de vorm: I = {m, m + 1, . . . , M },

A = {0, 1, 2, . . . , M − m}.

De oorspronkelijke toestandsruimte en actieruimte waren aftelbaar oneindig, maar aangezien we in een positief recurrent voorraadmodel zeer grote voorraden en zeer grote tekorten nauwelijks tegen zullen komen, kunnen we het ons veroorloven deze toestandsruimte af te kappen, mits M groot genoeg en m klein genoeg gekozen worden. Het stuk code dat hiervoor geschreven is, staat zeer grote M en kleine m toe, waardoor de resultaten zeer nauwkeurig zijn. In de volgende sectie geven we een oplossing van ons voorraadprobleem middels het SA-algoritme, toegepast op het Markov beslissingsproces dat we hebben opgebouwd in voorgaande secties.

4.5. Optimale bestelstrategie middels SA Middels MATLAB heeft de auteur een programma (zie Appendix (A)) geschreven dat het SA-algoritme toepast op dit specifieke voorraadprobleem. Dit programma definieert een functie SAvoorraad met de volgende argumenten, in volgorde: • klanten per tijdseenheid λ, • levertijd L, • voorraadkosten per eenheid in voorraad per tijdseenheid r, • vaste kosten per bestelling K, • boetekosten per eenheid tekort per tijdseenheid b, • gebruikte truncatie trunc van de toestandsruimte, • precisie van het SA-algoritme. Het argument trunc vereist wat opheldering. Deze variabele geeft het maximum M aan waarop we de toestandsruimte afkappen. In het programma wordt de ondergrens m dan gegeven door m = −M . Het programma past vervolgens het SA-algoritme toe en geeft de strategie f als vector na voldoende iteraties zodanig dat Mn − mn < . Het programma levert verhelderende resultaten. Bij de formulatie van ons model hebben we de hypothese opgesteld dat de optimale bestelstrategie een (s, S)-strategie zou zijn, waarbij een hoeveelheid S −i besteld werd wanneer de voorraad i lager was dan s.We voeren het programma uit voor verschillende combinaties van de parameters λ, L, r, K en b. Telkens kiezen we trunc = 20 (dus de toestandsruimte is {−20, −19, . . . , 19, 20}) 1 . en = 1000

43

Eerst bekijken we de combinatie λ = 2, L = 1, r = 1, K = 1 en b = 3, die we veel gebruikt hebben in onze voorbeelden. Het invullen van deze parameters levert de volgende strategie f op: f (−20) = f (−19) = 24, f (−18) = 23, f (−17) = f (−16) = 22, f (i) = f (i − 1) − 1, i = −15, −14, . . . , 0, f (i) = 0, i = 1, 2, . . . , 20. Deze strategie komt bijna volledig overeen met een (s, S)-strategie, waarbij s = 0 en S = 6. Alleen bij de laagste voorraadniveaus (-20 tot en met -17) wijkt de strategie hier iets van af. Verder zien we dat de het bestelpunt s = 0 zeer laag is. Waarschijnlijk zijn de boetekosten hier dus te laag ingesteld: zo laag, dat de strategie geplande tekorten heeft. Voeren we de boetekosten op naar b = 10, dan verkrijgen we de volgende strategie f : f (−20) = 27, f (−19) = f (−18) = 26, f (−17) = 25, f (−16) = 24, f (−15) = f (−14) = 23, f (i) = f (i − 1) − 1, i = −13, −12, . . . , 3, f (i) = 0, i = 4, 5, . . . , 20. Nu hebben we een strategie die grotendeels overeenkomt met een (s, S)-strategie met s = 3 en S = 9. Hier zorgen tekorten voor dusdanig hoge kosten dat er een bestelpunt relatief ver boven 0 is gekozen (merk op dat het verwachte klanten tijdens de levertijd slechts λL = 2 is. Weer zien we afwijkingen van deze strategie bij de lagere voorraadniveaus. Dit valt te wijten aan het feit dat voor de ondergrens van de nettovoorraad de overgangskansen zijn aangepast om de toestandsruimte te trunceren. Dit heeft als gevolg dat de boetekosten bij voorraden dicht bij de ondergrens worden onderschat (het tekort kan niet verder oplopen dan de ondergrens). Dit effect wordt sterker naarmate λL stijgt ten opzichte van de grootte van de toestandsruimte, zoals de volgende combinatie van parameters laat zien. Hier stellen we λ = 4 en b = 3. De rest van de parameters blijft gelijk. De berekende strategie is als volgt: f (−20) = f (−19) = f (−18) = 27, f (−17) = f (−16) = 26, f (−15) = 25, f (−14) = f (−13) = 24, f (i) = f (i − 1) − 1, i = −12, −11, . . . , 2, f (i) = 0, i = 3, 4, . . . , 20.

44

Hier wijken meer waarden van f (i) af van een (s, S)-strategie. In dit voorbeeld is de gekozen truncatie slechts vijf maal het verwachte aantal klanten in de levertijd. Hierdoor worden toestanden weggelaten die nog significant vaak voorkomen in het voorraadproces. Daarom is het verstandig de toestandsruimte uit te breiden. De volgende strategie is berekend met dezelfde parameters als het vorige voorbeeld, met de aanpassing dat trunc = 50. Dit levert de volgende strategie op: f (−50) = 58, f (−49) = f (−48) = 57, f (−47) = f (−46) = 56, f (−45) = 55, f (−44) = f (−43) = 54, f (−42) = 53, f (−41) = 52, f (−40) = 51, f (−39) = f (−38) = 50, f (i) = f (i − 1) − 1, i = −37, −36, . . . , 2, f (i) = 0, i = 3, 4, . . . , 50. Het effect blijft bestaan voor de eerste elementen van de functie, maar voor grotere nettovoorraden is de optimale strategie weer een (s, S)-strategie op plekken waar de strategie voor een lagere truncatie dat niet was. Dat de strategie iets afwijkt van een (s, S)-strategie voor lage toestanden maakt minder uit, aangezien deze toestanden nauwelijks voor zullen komen. Bovendien hebben we hier een strategie die 10 bestelt op het bestelpunt in plaats van 9. Het is voor nauwkeurige resultaten dus aan te raden de toestandsruimte groot te houden. We kunnen wel concluderen dat de hypothese juist is en dat de optimale strategie gegeven wordt door een (s, S)-strategie. Niet alleen dat, maar we hebben een programma gebouwd dat ons effectief het optimale bestelpunt s en referentieniveau S geeft.

45

5. Discussie In deze scriptie hebben we een analyse gemaakt van een voorraadmodel dat gemodelleerd is als Markovketen. We hebben het model in eerste instantie analytisch bekeken, door middel van technieken die bij Markovketens veelvuldig worden gebruikt. Als toevoeging hierop hebben we de matrix-geometrische methode toegepast, om de betreffende Markovketen, dat aftelbaar oneindig veel toestanden bevat, te kunnen beschouwen als een eindig-dimensionaal matrixprobleem. Er bleek geen gesloten formule te bestaan als oplossing van het matrixprobleem, dus we moesten ons wenden tot het benaderen van de oplossing. Dit maakte een analytische oplossing van het optimale bestelpunt en optimale bestelgrootte niet meer mogelijk, maar we waren wel in staat de gemiddelde kosten op lange termijn weer te geven voor verschillende combinaties van het bestelpunt en de bestelgrootte. Dit gaf ons inzicht in het verloop van de kosten als functie van het bestelpunt en de bestelgrootte. De beste resultaten werden behaald door ons model te formuleren als Markov beslissingsproces. Binnen de Markov beslisslingstheorie bestonden zeer effectieve algoritmes. Een van deze algoritmes werd ge¨ımplementeerd in een stuk MATLAB-code, zodat numeriek de beste bestelstrategie benaderd kon worden. Dit optimum had vervolgens een vorm die bekend staat als een (s, S)-strategie, zoals we verwacht hadden bij het formuleren van het model. Oorspronkelijk was het plan het model te beschouwen met constante levertijden. Er is een uitstap gemaakt naar exponentieel verdeelde levertijden omdat de voorraad zich met deze aanname gedraagt als een continue-tijds Markovketen, wat wiskundig beter beheersbaar is. Constante levertijden zijn realistischer in dit model dan exponentieel verdeelde levertijden, aangezien er in de praktijk geen grote variatie bestaat in levertijden. Uiteindelijk bleek het model met exponentieel verdeelde levertijden echter bevredigende resultaten te leveren, waardoor ervoor gekozen is deze aanname te houden. Hierna zouden we een model bekijken waarbij de levertijd een som is van meerdere exponentieel verdeelde stochasten. Vanwege de sterke wet van de grote aantallen zouden we hiermee het model met constante levertijden kunnen benaderen, en misschien zelfs een beter model verkrijgen dat enige variatie in de levertijd toelaat. Aan dit proces zijn we echter niet toegekomen, gezien de tijd en inhoud dat de analyse van het model met exponentieel verdeelde levertijden in beslag nam. Verder onderzoek zou gedaan kunnen worden in dit vlak. Waar ten slotte nog aandacht aan kan worden besteed in vervolgonderzoek is het Periodic Review model. Dit model heeft veel weg van het gebruikte model, in die zin dat niet op elk tijdstip een bestelling geplaatst kan worden. De optimale strategie in zo’n model is eveneens een (s, S)-strategie. De overeenkomsten en verschillen van deze twee modellen zouden onderzocht kunnen worden.

46

Bibliografie [1] H.C. Tijms, E.M.F. Kalvelagen, Modelbouw in de operations research, Academic Service, Schoonhoven, 1994. [2] R. N´ un ˜ez-Queija, Markov Decision Processes - Introduction, Amsterdam, 2011. [3] H.C. Tijms, John Wiley & Sons, Stochastic Models: An Algorithmic Approach, Chichester, 1994. [4] Prez, J. F.; Van Houdt, B. (2011), Quasi-birth-and-death processes with restricted transitions and its applications, Performance Evaluation 68 (2). [5] M.F. Neuts, Matrix-Geometric Solutions in Stochastic Models: An Algorithmic Approach, The John Hopkins University Press, Baltimore, 1981.

47

A. SA-algoritme in MATLAB function f=SAvoorraad(lambda,L,r,K,b,trunc,epsilon) I=trunc*2+1; A=I; X=I; R=zeros(A,I); toestanden=-trunc:trunc; gamma=L/(lambda*L+1); for i=1:trunc R(1,i)=-gamma*toestanden(i)*b; end for i=trunc+1:I R(1,i)=gamma*toestanden(i)*r; end for i=1:trunc R(2:A,i)=(-gamma*toestanden(i)*b+K)*ones(A-1,1); end for i=trunc+1:I R(2:A,i)=(gamma*toestanden(i)*r+K)*ones(A-1,1); end v=zeros(I,X,2); f=zeros(1,I); for i=1:I v(i,1,2)=R(1,i); f(i)=1; if i
48

M=max(v(1,:,2)); m=min(v(1,:,2)); for i=2:I if max(v(i,I+1-i,2))>M M=max(v(i,I+1-i,2)); end if min(v(i,I+1-i,2))<m m=min(v(i,I+1-i,2)); end end pk=lambda*L/(lambda*L+1); pl=1/(lambda*L+1); while M-m>=epsilon v(:,:,1)=v(:,:,2); v(1,1,2)=R(1,1)+v(1,1,1); f(1)=1; for a=2:A if R(a,1)+pk*v(1,a,1)+pl*v(a,1,1)
49

f(I)=1; M=max(v(1,:,2)-v(1,:,1)); m=min(v(1,:,2)-v(1,:,1)); for i=2:I if max(v(i,I+1-i,2)-v(i,I+1-i,1))>M M=max(v(i,I+1-i,2)-v(i,I+1-i,1)); end if min(v(i,I+1-i,2)-v(i,I+1-i,1))<m m=min(v(i,I+1-i,2)-v(i,I+1-i,1)); end end end f=f-ones(1,I);

50

B. Populaire samenvatting Managen van een bedrijf is keuzes maken. Keuzes omtrent het personeel, omtrent marketing, omtrent prijszetting. Er zijn vele motiveringen te vinden voor een keuze. We hebben hier te maken met mensen, dus de integriteit van het bedrijf is een belangrijke factor. Maar voor het gros van de keuzes zijn beweegredenen te reduceren tot de volgende centrale vraag: Hoe maximaliseer ik mijn winst? De voor de hand liggende oplossing van deze vraag is het maximaliseren van opbrengsten. Winst wordt echter niet alleen bepaald door de opbrengsten die een bedrijf heeft, maar ook door de kosten. De voorraadtheorie houdt zich bezig met een specifieke vorm van kosten. De naam zegt het al: het gaat hier om het managen van de voorraad van een product. De kosten die een bedrijf maakt bij het houden van voorraad kunnen worden gezien in de volgende schets. Stel je bent manager van het bedrijf TV Discount Online. Je verkoopt vele soorten TV’s, maar de nieuwste 3D-televisie van Samsung is razend populair geworden. Daarom besluit je een voorraad van deze TV’s te houden. De TV’s zijn zeer waardevol, dus je besluit verzekeringen af te sluiten op de TV’s die je inkoopt, voor het geval dat er een brand uitbreekt in het magazijn die alle televisies verwoest. Je weet hoeveel klanten gemiddeld per dag de TV bestellen. Soms zijn het er meer, soms minder. Omdat de TV zo populair is, willen klanten hem zo snel mogelijk in huis hebben. Om klanten niet aan een concurrent te verliezen als de voorraad op is, geef je een compensatie aan elke klant die groter wordt hoe langer ze op hun product moeten wachten. Uiteraard moet er worden bijbesteld als de voorraad leeg dreigt te raken. Je doet dan een bestelling bij de leverancier, die een vrachtwagen stuurt met het aangegeven aantal TV’s. De vervoersen administratiekosten moet je zelf betalen. De manager in dit voorbeeld heeft te maken met verschillende soorten kosten. Ten eerste heeft hij kosten aan het in voorraad houden van de TV’s (de verzekeringen). Deze kosten stijgen hoe meer hij in voorraad houdt. Ten tweede zijn er vaste kosten (vervoersen administratiekosten) voor elke bestelling. Hoe vaker hij bestelt, hoe hoger deze kosten worden. Ten slotte heeft hij kosten aan het uitverkocht zijn (de boetekosten). Deze kosten worden groter hoe meer en hoe langer klanten tegen een lege voorraad aanlopen. De manager heeft nu twee beslissingen om te maken: • Hoeveel moet je bestellen? • Wanneer moet je bestellen? Deze twee vragen zijn niet gemakkelijk te beantwoorden. Bestelt de manager te veel, dan komt hij te zitten met hoge voorraadniveaus en dus hoge voorraadkosten. Bestelt hij te weinig, dan moeten er te veel kleine bestellingen worden gedaan, waardoor er te veel betaald wordt aan vaste kosten. Bestelt de manager te laat, dan raakt de voorraad op

51

terwijl de vrachtwagen onderweg is en moeten klanten wachten op hun product. Bestelt de manager te vroeg, dan heeft hij weer hoge voorraadkosten. Ergens ligt een optimum voor de hoeveelheid en het tijdstip van bestellen, dat de kosten minimaliseert. Het is de taak aan de voorraadtheorie dit optimum te vinden. Er zijn verscheidene modellen die dit probleem behandelen. In deze scripte behandelen we een model uitvoerig. We kijken hoe de keuze van de bestelgrootte en tijdstip van bestellen de kosten op lange termijn be¨ınvloedt. De voorraad wordt gemodelleerd als Markovketen, wat wil zeggen dat de voorraad op een willekeurige manier in de tijd verandert. De scriptie sluit af met de oplossing van het probleem, berekend door de computer. Deze oplossing blijkt een zeer bekende vorm te hebben, die veel gebruikt wordt in andere modellen.

52

Voorraadtheorie. Mathijs van der Vlies 28 augustus Begeleiding: prof.dr. R. Núñez Queija

Recommend Documents