Even geduld a.u.b. 3 maart 2003 In het dagelijks leven hebben we vaak te maken met wachten. Denk bijvoorbeeld aan het wachten bij de kassa in de supermarkt, voor het downloaden van een file op het internet en voor stoplichten. Wachten vinden we meestal niet prettig. Om het wachten te verkorten zijn aanpassingen nodig (zoals extra kassa’s), maar de vraag is: Wat en hoe? Om te besluiten wat nodig is willen we graag antwoord op vragen zoals: Wat is de wachttijd van een klant als we een of meer extra kassa’s openen? Hoe lang duurt het om een file te downloaden als we een snelle internetverbinding hebben? Een antwoord op deze vragen kan worden gegeven met behulp van wiskundige modellen die wachtrijsituaties beschrijven: wachtrijmodellen. De oorsprong van de wachtrijtheorie gaat terug naar het begin van de vorige eeuw toen de Deense ingenieur A.K. Erlang (1878-1929) een studie publiceerde over de kansverdeling van het aantal bezette telefoonlijnen in een bundel van N lijnen (zie figuur 1 en de referenties [2, 3]). Dit werk is nog steeds relevant: de Erlang-B formule wordt bijvoorbeeld toegepast bij het ontwerp van mobiele telecommunicatienetwerken (van bijvoorbeeld KPN en Vodafone), en is ook bruikbaar bij het ontwerp van parkeergarages, databases en containeroverslagterreinen (zoals van ECT in Rotterdam). Toepassingen van wachtrijmodellen voor het ontwerp van computersystemen, communicatienetwerken, produktiesystemen en call centers worden bijvoorbeeld beschreven in de referenties [6, 11, 1, 4, 5]; een minder gebruikelijke toepassing is de applet Cow waarin een wachtrijmodel wordt gebruikt voor het ontwerp van een moderne koeienstal met melkrobots. Bij het bestuderen van wachtrijmodellen speelt de kansrekening een belangrijke rol, want wachten wordt juist veroorzaakt door toevallige schommelingen in aankomsttijden en bedieningstijden. Maar kansrekening is een lastig onderwerp om je eigen te maken. De computer biedt nieuwe mogelijkheden om ons te helpen wachtrijsituaties te begrijpen. Bij computersimulatie wordt een wachtrijsituatie nagespeeld op de computer en kunnen resultaten op heldere wijze op het scherm worden getoond. Een uitstekende inleiding over kansrekening en computersimulatie wordt gegeven in het boek [8], en ook in het veel uitgebreidere boek [10]. In deze Masterclass gaan we kennismaken met wachtrijmodellen en computersimulatie. Dit zal gebeuren aan de hand van twee herkenbare praktijksituaties: het bepalen van een stoplichtafstelling en het ontwerpen van een parkeerterrein voor een grote supermarkt. Maar 1
1 N ρ N! B(N ) = 1 1 N 1 + ρ + ρ2 + · · · + ρ 2 N!
Figuur 1: De Deense ingenieur A.K. Erlang (1878-1929) en zijn formule misschien wel belangrijker is dat we naast wachtrijtheorie en simulatie kennis zullen zaken met Modelleren. Dit is het vertalen van een praktijkprobleem naar een wiskundig probleem. Vragen zoals ‘Hoe vinden we een betere stoplichtafstelling?’ of ‘Hoe groot dient het parkeerterrein te zijn?’ zijn helder, maar niet wiskundig. Het bedenken van een passend wiskundig probleem zullen we zelf moeten doen; wiskundige problemen komen nu eenmaal niet in het wild voor! Het zal duidelijk zijn dat modelleren zeer belangrijk is voor een wiskundig ingenieur. In het eerste deel van de masterclass zullen we ons richten op het stoplichtprobleem; in het tweede deel komt de parkeerplaats aan de orde.
1
Afstelling van stoplichten
We hebben allemaal de neiging om te mopperen op de afstelling van stoplichten: dat kan toch veel beter! Maar het vinden van een goede afstelling van stoplichten is helemaal niet eenvoudig. Om enig gevoel voor dit probleem te krijgen zullen we gaan kijken naar de volgende situatie met slechts twee verkeersstromen. Wegwerkzaamheden Wegens wegwerkzaamheden is een van de twee rijstroken van een tweebaansweg over een afstand van 500 meter afgesloten (zie figuur 2). Het verkeer uit beide richtingen wordt door middel van stoplichten om beurten over de vrije rijstrook geleid. Er geldt een snelheidsbeperking van 20 kilometer per uur. De vraag is om een goede afstelling voor de stoplichten te bepalen. Om een goede afstelling te bepalen helpt het enorm als we de verkeersituatie kunnen vertalen naar een wiskundig model. Dit model dient de situatie voldoende nauwkeurig te 2
Figuur 2: Wegwerkzaamheden over 500 m beschrijven en kan dan worden gebruikt om wachttijden en rijlengtes te voorspellen voor een gegeven afstelling. De voorspellingen kunnen worden gebruikt voor het vinden van een ‘optimale’ afstelling: in dit geval ligt het voor de hand om gemiddelde wachttijd te minimaliseren. Het is natuurlijk de kunst om een goed model op te stellen. In dit voorbeeld kijken we naar gewoon verkeer, maar op het internet treffen we een soortgelijke situatie aan: op het internet zijn vele computers via een (verkeers)netwerk met elkaar verbonden (zie figuur 3). Het netwerk bestaat uit links (de verkeerswegen) die met elkaar zijn verbonden door zogenaamde switches (de kruispunten) en hierover versturen computers stromen datapakketjes of files (het verkeer) naar elkaar.
Figuur 3: Computercommunicatienetwerk
1.1
Modellering
Laten we de stroom auto’s in ´e´en richting bekijken. De stroom in de andere richting kan op precies dezelfde manier behandeld worden. Voor een stoplichtafstelling is het van belang hoeveel auto’s arriveren en hoe snel ze het stoplicht passeren. De aankomstintensiteit van auto’s geven we aan met λ auto’s per tijdseenheid en de vertrekintensiteit van auto’s 3
Figuur 4: Verloop van het aantal wachtende auto’s N (t) als functie van de tijd voor λ = 1/3, µ = 1, r = 9 en g = 6 gedurende groentijden met µ auto’s per tijdseenheid; dus 1/λ is de gemiddelde tijd tussen twee (opeenvolgende) aankomsten en 1/µ de gemiddelde tijd tussen twee passerende auto’s. We nemen aan dat zowel de tussenaankomst- als passeertijden constant zijn (hier komen we later nog op terug). Als tijdseenheid kunnen we bijvoorbeeld nemen uur, minuut of seconde. In het vervolg kiezen we seconde als tijdseenheid. Verder is r de roodtijd (in seconden) en g de groentijd in een cyclus; de cyclustijd geven we aan met c, dus c = r + g. Natuurlijk moet gelden λ · c ≤ µ · g. (1) Aan de linkerkant staat namelijk het aantal auto’s dat in een cyclus aankomt, terwijl aan de rechterkant het aantal auto’s staat dat maximaal in een cyclus kan vertrekken; dus als niet aan bovenstaande voorwaarde is voldaan, zal het aantal wachtende auto’s onbeperkt groeien. Onder de voorwaarde (1) kunnen we de gemiddelde wachttijd van auto’s bij het stoplicht uitrekenen. Allereerst zullen we kijken naar het gemiddeld aantal auto’s dat bij het stoplicht staat te wachten. Gedurende een cyclus varieert het aantal auto’s dat staat te wachten op de volgende wijze (zie figuur 4): Het aantal wachtende auto’s verloopt als een stapfunctie. Voor het gemak benaderen we de stapfunctie door een continue functie; gedurende de roodtijd stijgt deze functie met richtingsco¨efficient λ en gedurende de groentijd daalt deze functie met richtingsco¨efficient λ − µ. De oppervlakte O onder de (continue) grafiek is gelijk aan O=
1 λr λr2 1 · r · λr + · · λr = , 2 2 µ−λ 2(1 − ρ)
waarin ρ = λ/µ. Het gemiddeld aantal auto’s N dat staat te wachten wordt dan gegeven door de oppervlakte onder de grafiek gedeeld door de cycluslengte, dus N=
λr2 O = c 2c(1 − ρ)
(2)
Om nu de gemiddelde wachttijd W van auto’s uit te rekenen maken we gebruik van een beroemde formule uit de wachtrijtheorie, namelijk de formule van Little (zie referentie [7]). Deze formule luidt N = d · W, (3) 4
waarin • N het gemiddeld aantal wachtende klanten in het wachtrijsysteem is; • d de doorzet van het wachtrijsysteem is, dat wil zeggen het gemiddeld aantal klanten dat het systeem per tijdseenheid afhandelt; • W de gemiddelde wachttijd van een klant in het wachtrijsysteem is. De volgende intu¨ıtieve redenering verklaart waarom deze formule plausibel is. Stel dat klanten voor de tijd dat ze wachten in het wachtrijsysteem 1 euro per tijdseenheid krijgen van de beheerder van het systeem. Er zijn dan twee mogelijkheden waarop de beheerder dit geld kan uitbetalen: • Iedere klant krijgt bij vertrek het totale bedrag uitbetaald. Met een doorzet van d klanten per tijdseenheid en een gemiddelde wachttijd van W per klant, zal de beheerder gemiddeld d · W euro per tijdseenheid uitbetalen. • De beheerder betaalt iedere tijdseenheid alle op dat moment wachtende klanten 1 euro. Aangezien er gemiddeld N wachtende klanten zijn, zal hij dus ook gemiddeld N euro per tijdseenheid uitbetalen. De beheerder zal in beide gevallen evenveel uitbetalen en dus geldt (3). We kunnen concluderen dat de gemiddelde wachttijd van auto’s bij ons stoplicht gegeven wordt door r2 (c − g)2 N = = . (4) W = λ 2c(1 − ρ) 2c(1 − ρ)
Voorbeeld 1.1 Voor het voorbeeld uit figuur 4 vinden we (ga na) W =
1.2
81 (15 − 6)2 = (seconden). 2 · 15 · (1 − 1/3) 20
(5)
De optimale afstelling
Voor het vinden van de optimale afstelling dienen we rekening te houden met beide verkeersstromen. De aankomstintensiteit van auto’s van stroom i geven we aan met λi en de groen- en roodtijd voor stroom i met gi en ri (i = 1, 2). Nieuw is de ontruimingstijd van de weg, aangegeven met x. Dit is de minimale tijd die verstrijkt vanaf het moment dat het ene stoplicht op rood springt totdat het andere stoplicht op groen springt. Met andere woorden, x is (tenminste) de tijd die nodig is om 500 meter af te leggen bij een snelheid van 20 km/uur, en gedurende deze tijd staan beide stoplichten op rood. De cyclustijd c van het stoplicht kunnen we schrijven als c = g1 + r1 = g2 + r2 = g1 + x + g2 + x. 5
Uiteraard moet weer gelden (zie (1)) λi · c ≤ µ · gi ,
i = 1, 2.
Onder deze voorwaarden wordt de gemiddelde wachttijd Wi van de auto’s van stroom i gegeven door (zie (4)) (c − gi )2 Wi = , i = 1, 2, (6) 2c(1 − ρi ) waarin ρi = λi /µ. Voor de gemiddelde wachttijd W (g1 , g2 ) van alle auto’s bij groentijden g1 en g2 vinden we λ2 λ1 W1 + W2 . (7) W (g1 , g2 ) = λ 1 + λ2 λ1 + λ 2 Het vinden van de optimale stoplichtafstelling kan nu worden vertaald in het volgende (wiskundige) optimaliseringsprobleem: Minimaliseer de functie W (g1 , g2 ) onder de voorwaarden (i) gi ≥ 0, (ii) λi · c ≤ µ · gi ,
i = 1, 2; i = 1, 2.
Dit lijkt wat op een lineair programmeringsprobleem. De voorwaarden zijn namelijk lineair; echter, de doelfunctie is dat niet. De voorwaarden onder (ii) kunnen ook herschreven worden tot λ2 λ2 g1 + 2x, µ − λ2 µ − λ2 µ − λ1 ≤ g1 − 2x. λ1
g2 ≥ g2
Het toegelaten gebied is dus van de volgende vorm: Merk op dat de kleinst mogelijke waarden van g1 en g2 die we kunnen kiezen gegeven worden door het snijpunt van de twee lijnen in bovenstaande figuur, dus λ1 · 2x, µ − λ1 − λ 2 λ2 = · 2x. µ − λ1 − λ 2
g1 = g2
De bijbehorende (kleinste) cyclustijd is c=
µ · 2x. µ − λ 1 − λ2
Maar zijn dit ook de best mogelijke waarden van g1 en g2 die we kunnen kiezen? Om dat te kunnen bepalen gaan we het optimaliseringsprobleem oplossen in twee stappen: 6
−
−
−
−
Figuur 5: Toegelaten gebied voor de groentijden g1 en g2 (1) Allereerst bepalen we voor elke vaste waarde van de cycluslengte c de best mogelijke keus van g1 en g2 ; (2) Daarna bepalen we de best mogelijke keus van c. Hoewel er veel meer te zeggen is over waar de optimale waarden van g1 en g2 liggen, doen we dat niet, maar gaan we bovenstaande stappen op de computer uitvoeren met behulp van de volgende applet StoplichtSim.
Voorbeeld 1.2 Stel λ1 = 1/3, λ2 = 1/4, µ = 1 en x = 2. Dus ρ1 = 1/3, ρ2 = 1/4, c = g1 + g2 + 4 en W1 =
3(g2 + 4)2 , 4(g1 + g2 + 4)
W2 =
3(g1 + 4)2 . 2(g1 + g2 + 4)
Het minimaliseringsprobleem kan nu worden geschreven als Minimaliseer de functie W (g1 , g2 ) =
4 3(g2 + 4)2 3 3(g1 + 4)2 · + · 7 4(g1 + g2 + 4) 7 2(g1 + g2 + 4)
onder de voorwaarden (i) g1 ≥ 0, (ii) g2 ≥ 21 g1 + 2,
g2 ≥ 0, g2 ≤ 2g1 − 4.
De kleinst mogelijke groentijden zijn g1 = 3 51 en g2 = 2 25 .
7
Voor het speciale geval dat beide verkeersstromen dezelfde intensiteit hebben kunnen we de optimale groentijden direct (zonder computer) bepalen. Dit is het doel van de eerste opdracht. Opgave 1. Neem aan dat de intensiteiten van beide stromen hetzelfde zijn, dus λ1 = λ2 = λ. Het ligt voor de hand dat ook de optimale groen- en roodtijden voor beide richtingen hetzelfde zijn, dus g1 = g2 = g en r1 = r2 = r. (i) Laat zien dat de gemiddelde wachttijd W (g) voor alle auto’s gelijk is aan (c − g)2 W (g) = , 2c(1 − ρ) waarin c = 2(g + x) en ρ = λ/µ. Neem nu aan dat λ = 1/3, µ = 1 en x = 2. (ii) Ga na dat W (g) =
3(g + 4)2 . 8(g + 2)
(iii) Toon aan dat W (g) stijgend is voor g ≥ 0. (iv) Laat zien dat uit (ii) volgt dat de optimale groentijd gelijk is aan de minimale groentijd g = 4. In grafiek 6 staan de gemeten aankomstintensiteiten voor de verkeersstromen uit beide richtingen. Zoals je kunt zien is er sprake van een ochtend- en avondspits. Bij de volgende opdracht gebruiken we de applet Naam voor het bepalen van de optimale groentijden tijdens de spits en daluren. Een redelijke waarde voor de passeertijd is 2 seconden; ga maar eens meten bij een stoplicht! Opgave 2. (i) Welke spits- en daluren wil je onderscheiden op basis van grafiek 6? En wat zijn gedurende deze uren de aankomstintensiteiten? (ii) Bereken de ontruimingstijd en de minimale cyclusduren gedurende de spits- en daluren. (iii) Bepaal voor zowel de spits- als daluren de optimale groentijden met behulp van de applet Naam. (vi) Probeer de verschillen in de cycli gevonden onder (iii) te verklaren.
8
1000
800
600
400
200
0 6
8
10
12
14
16
18
20
Figuur 6: Gemeten aankomstintensiteiten (auto’s per uur) voor stroom 1 (doorgetrokken lijn) en stroom 2 (onderbroken lijn) van 6.00 uur tot 20.00 uur
1.3
Fluctuaties in aankomsten
Laten we het model nog eens onder de loep nemen. We hebben aangenomen dat zowel de passeertijden als aankomsttijden constant zijn. Dit lijkt redelijk voor de passeertijden, maar is veel minder realistisch voor de aankomsttijden van auto’s. Naast (geleidelijke) veranderingen in de aankomstintensiteit kunnen in de praktijk de tussenaankomsttijden van auto’s flink schommelen. Een model dat (vaak) een goede beschrijving geeft van de schommelingen in de tussenaankomsttijden is het Poisson proces (zie bijvoorbeeld de referentie [9]). In dit model veronderstellen we dat de tijd die verstrijkt tussen twee opeenvolgende aankomsten exponentieel verdeeld is en onafhankelijk is van alle andere tussenaankomsttijden. Als we een tussenaankomsttijd aangeven met A (we noemen dit een stochastische variabele), dan nemen we dus aan dat F (t) = P (A ≤ t) = 1 − e−λt ,
t ≥ 0.
We noemen dit ook de exponenti¨ele verdeling. In figuur 7 is de verdelingsfunctie F (t) weergegeven; we zien bijvoorbeeld dat F (1.3) = 0.73, ofwel 73% van de tussenaankomsttijden is kleiner of gelijk aan 1.3 seconden. Figuur 8 laat zien hoe sterk de tussenaankomsttijden van een Poissonproces kunnen schommelen. Een Poissonproces betekent in feite dat de aankomsten volkomen toevallig zijn. Het model waarin auto’s arriveren volgens een Poisson proces zal de verkeerssituatie beter beschrijven. Een nadeel is dat het veel moeilijker is om te analyseren. Maar het is natuurlijk denkbaar dat ons eenvoudige model in sectie 1.1 een prima stoplichtafstelling levert. We kunnen dit direct onderzoeken met behulp van een computersimulatie. Dat wil zeggen, we 9
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
1
2
3
4
Figuur 7: De verdelingsfunctie F (t) met λ = 1
Figuur 8: Aankomsten volgens een Poissonproces
10
5
gaan de verkeerssituatie naspelen op de computer. Dan kunnen we bekijken of de wachttijden die de computer produceert overeenkomen met onze schattingen gevonden in Opgave 2. In het computermodel moeten we schommelende tussenaankomsttijden genereren. Preciezer gezegd: we dienen random trekkingen te genereren uit een exponenti¨ele verdeling. Hieronder beschrijven we hoe je dat op een computer kunt realiseren. Op elke computer is tegenwoordig een random generator ingebouwd. Een random generator kiest getallen uit die willekeurig verspreid liggen tussen 0 en 1. Zo’n getal heet een random getal of ook wel toevalsgetal. Een random generator volstaat om trekkingen te doen uit vrijwel elke denkbare kansverdeling. In de volgende opdracht is het de bedoeling dat jullie dit laten zien voor de exponenti¨ele verdeling. Opgave 3. Zij U een willekeurig gekozen getal tussen 0 en 1, dus P (U ≤ x) = x,
0 ≤ x ≤ 1.
Bekijk nu de stochastische variabele A = F −1 (U ) =
− ln(1 − U ) . λ
Toon aan dat P (A ≤ t) = 1 − e−λt ,
t ≥ 0,
ofwel, A is exponentieel verdeeld! Dus met de computer kunnen we als volgt een random getal a trekken uit de exponenti¨ele verdeling; random is de random generator (procedure) van de computer en lambda is de parameter van de exponentiele verdeling. 1. u := random 2. a := - ln(1 - u) / lambda 3. return a In grafiek 7 is bijvoorbeeld te zien dat als de computer het getal u = 0.73 kiest, dan is de bijbehorende tussenaankomsttijd a = 1.3. Het bovenstaand algoritme voor een trekking uit de exponenti¨ele verdeling is eenvoudig te implementeren op je grafische rekenmachine (ga na). Nu gaan we de stoplichtafstellingen gevonden in Opdracht 2 simuleren met behulp van de applet Naam. In het simulatiemodel gebruiken we de (realistische) Poisson beschrijving voor aankomsten van auto’s. Opgave 4.
11
(i) De computersimulatie zal schattingen voor gemiddelde wachttijden en rijlengtes produceren die (een beetje) schommelen, want de auto’s komen op willekeurige momenten aan. Maar als je voldoende lang simuleert krijg je nauwkeurige schattingen. Wat is voldoende lang? (ii) Simuleer met behulp van de applet Naam de stoplichtafstellingen voor spits- en daluren gevonden in Opdracht 2. (ii) Beschrijf je bevindingen met computersimulatie. Leveren de afstellingen uit Opdracht 2 goede resultaten? (iii) Probeer met behulp van computersimulatie betere afstellingen te vinden.
1.4
De formule van Webster
Je kunt waarschijnlijk met de computersimulatie door zomaar proberen een goede afstelling voor de stoplichten vinden. Maar voor problemen met veel verkeersstromen (denk bijvoorbeeld aan een groot kruispunt) is dit al snel ondoenlijk. Het is veel effici¨enter om een goede afstelling te vinden op basis van een wiskundige aanpak. Het model met constante tussenaankomsttijden van auto’s is eenvoudig, maar geeft blijkbaar (te) slechte schattingen voor de wachttijden. Anderzijds geeft het model met variabele (exponenti¨ele) tussenaankomsttijden een betere beschrijving, maar het is veel lastiger. Voor dit model zullen we nu een formule geven waarmee de gemiddelde wachttijd van auto’s vrij nauwkeurig voorspeld kan worden. Deze formule heet de formule van Webster (zie referentie [12]), en met behulp van deze formule zullen we proberen om wel een goede afstelling te vinden. Laten we eerst de stroom auto’s in ´e´en richting bekijken (net zoals in sectie 1.1). We nemen aan dat de tussenaankomsttijden exponentieel verdeeld zijn en dat de passeertijden constant zijn. Voor de groentijd g in een cyclus moet natuurlijk weer gelden (zie (1)) µ · g ≥ λ · c. Een goede benadering voor de gemiddelde wachttijd W van auto’s kan worden verkregen met de formule van Webster, W =
ρc2 (c − g)2 + . 2c(1 − ρ) 2g(µg − λc)
(8)
Voor het voorbeeld uit figuur 4 levert dit (vergelijk (5) W =
81 25 + (seconden). 20 4
Als we formule (8) bekijken, zien we dat de eerste term overeenkomt met formule (4). Deze term is dus de gemiddelde wachttijd van auto’s als de tussenaankomsttijden constant 12
zouden zijn, en is in feite de bijdrage van de stoplichtafstelling. De tweede term is de bijdrage van schommelingen in de tussenaankomsttijden; deze term kunnen we als volgt uitleggen. We kunnen de roodtijd verwerken in de passeertijd door deze op te rekken met een factor c/g; de opgerekte passeertijd wordt 1 1 c = · . τ µ g Merk op dat voor opgerekte passeertijden het maximaal aantal auto’s dat het stoplicht kan passeren in een cyclus c gelijk is aan c · τ = µg. Dit is natuurlijk precies hetzelfde is als in de oorspronkelijke situatie. Een arriverende auto treft gemiddeld L wachtende auto’s aan bij het stoplicht plus (mogelijk) een voorste auto die bezig is om het stoplicht te passeren. De restpasseertijd van de voorste auto is gemiddeld de helft van een volledige passeertijd, en de kans om bij aankomst een (voorste) auto aan te treffen is gelijk aan λ · 1/τ . Na het passeren van de voorste auto duurt het nog L · 1/τ tijdseenheden totdat de wachtende auto’s het stoplicht zijn gepasseerd. We vinden dus W =
1 λ 1 λ 1 · +L· = 2 +L· . τ 2τ τ 2τ τ
Daarnaast geldt ook de formule van Little, dus L = λ · W. We hebben nu twee vergelijkingen voor twee onbekenden, W en L. Substitutie van bovenstaande formule in de formule voor W levert W =
λ/(2τ 2 ) ρc2 = , 1 − λ/τ 2g(µg − λc)
en dit is precies de tweede term in de formule van Webster! In grafiek 9 laten we de gemiddelde wachttijd volgens de formule van Webster zien als functie van de groentijd. We begrijpen nu waarom de eerste term uit de formule van Webster zo’n slechte schatting voor de gemiddelde wachttijd W geeft: als de groentijd g net iets groter is dan de minimale waarde λc/µ, ofwel µg − λc heel klein is, dan is de tweede term uit de formule van Webster veel groter dan de eerste term. Dit explosieve gedrag van de wachttijd als er maar net voldoende capaciteit is (in dit geval: groentijd), is typerend voor vrijwel elke wachtrijsituatie! We gaan nu verder met het vinden van de optimale stoplichtafstelling, waarbij we rekening houden met schommelde tussenaankomsttijden voor beide verkeersstromen. Met behulp van de formule van Webster vinden we voor de gemiddelde wachttijd W (g1 , g2 ) van alle auto’s (vergelijk formules (6) en (7)) W (g1 , g2 ) =
λ1 λ2 W1 + W2 λ1 + λ2 λ1 + λ2 13
20
15
10
5
0 5
6
7
8
9
10
11
12
13
14
15
Figuur 9: Gemiddelde wachtijd volgens Webster als functie van de groentijd voor λ = 1/3, µ = 1 en c = 15; de gestippelde lijn is de eerste term uit formule (8), de onderbroken lijn de tweede en de doorgetrokken lijn de som van de twee waarin
(c − gi )2 ρi c 2 + , i = 1, 2. 2c(1 − ρi ) 2gi (µi gi − λi c) Het bepalen van de optimale stoplichtafstelling kan weer worden vertaald in het optimalseringsprobleem uit sectie 1.2, maar dan met bovenstaande doelfunctie W (g1 , g2 ). Dit probleem gaan we weer oplossen met behulp van de computer. Wi =
Opgave 5. (i) Bepaal voor de spits- als daluren uit Opdracht 2 de optimale groentijden met behulp van de applet Naam. (ii) Simuleer de gevonden stoplichtafstellingen met de applet Naam en vergelijk de uitkomsten voor de wachttijden met de resultaten uit (i). (iii) Kun je de stoplichtafstellingen nog wat verbeteren? Opmerking 1.3 Eigenlijk is de formule van Webster iets ingewikkelder dan formule (8). Er hoort nog een derde (vreemde) term bij, namelijk c 1/3 λc 2+5g/c ρc2 (c − g)2 W = + − 0.65 2 . 2c(1 − ρ) 2g(µg − λc) λ µg 14
Deze derde term is een correctie term en is op een verstandige manier experimenteel bepaald, zodat de benadering voor de gemiddelde wachttijd nog iets beter is. De derde term is in het algemeen klein ten opzichte van de andere twee, en daarom hebben we deze term maar weggelaten.
2
Het ontwerpen van een parkeerterrein
Wanneer je wil gaan winkelen is het heel vervelend als je bij aankomst met de auto ontdekt dat er geen plaats is op het parkeerterrein bij het winkelcentrum. Maar hoe groot dient het parkeerterrein te zijn opdat dit nauwelijks voorkomt? We gaan nu proberen om op deze vraag een antwoord te geven. Parkeerterrein Een XL winkelcentrum heeft besloten om het parkeerterrein voor winkelende klanten te vergroten, zodat het plaats kan bieden voor vrijwel alle klanten. Tijdens drukke dagen (zoals vrij- en zaterdagen) verwacht men 200 klanten per uur, die gemiddeld 1.5 uur winkelen. Hoe groot dient het parkeerterrein te zijn? Om een antwoord op bovenstaande vraag te kunnen geven moeten we ons eerst afvragen wat het ontwerpcriterium is. We gaan er vanuit dat het percentage klanten dat bij aankomst een vol parkeerterrein aantreft voldoende klein dient te zijn; dit percentage noemen we het verliespercentage (teleurgestelde klanten besluiten wellicht om ergens anders te gaan winkelen). Wat we hier bedoelen met voldoende klein is iets dat uiteraard afgestemd moet worden met de leiding van het winkelcentrum. In de volgende sectie gaan we de situatie op het parkeerterrein wiskundig modelleren en met behulp van dit model kunnen we dan, afhankelijk van de grootte van het parkeerterrein, het verliespercentage voorspellen. We richten ons hier op de vraag of voldoende klanten een parkeerplaats kunnen vinden, maar een soortgelijke vraag kom je bijvoorbeeld tegen bij het ontwerp van een mobiel communicatienetwerk. Als je een telefoongesprek met je mobieltje wil voeren, dan dient de dichtsbijzijnde zendmast (het parkeerterrein) hiervoor een frequentie (een parkeerplaats) waarover het gesprek zal worden gevoerd vrij te hebben; zoniet, dan komt de verbinding niet tot stand of er wordt gekeken of nabij gelegen zendmasten wel frequenties beschikbaar hebben. Een belangrijke vraag is dus hoeveel (en welke) frequenties een zendmast dient te hebben opdat vrijwel alle gesprekken tot stand kunnen worden gebracht.
2.1
Modellering
In het model voor het parkeerterrein moeten we beschrijven: (i) wanneer klanten arriveren, (ii) hoelang ze gaan winkelen, en 15
Figuur 10: Fluctuaties in het aantal auto’s op het parkeerterrein (met capaciteit Max) (iii) wat ze doen als bij aankomst het parkeerterrein vol blijkt te zijn. Zowel de tussenaankomsttijden als de winkeltijden van klanten zullen (flinke) schommelingen vertonen; klanten kunnen op willekeurige momenten aankomen en ze gaan zeker niet allemaal even lang winkelen. De aankomsten van klanten (of auto’s) kunnen we net zoals voor het stoplichtprobleem beschrijven met een Poissonproces; de tussenaankomsttijden zijn exponentieel verdeeld met intensiteit λ klanten per uur. Verder veronderstellen we dat ook de winkeltijden van klanten exponentieel verdeeld zijn. De gemiddelde winkeltijd is 1/µ uur. Als we een tussenaankomsttijd aangeven met A en een winkeltijd met B, dan nemen we dus aan dat P (A ≤ t) = 1 − e−λt ,
P (B ≤ t) = 1 − e−µt ,
t ≥ 0.
Tot slot gaan we er vanuit dat klanten die bij aankomst een vol parkeerterrein aantreffen onmiddellijk vertrekken en niet meer terugkomen of ergens in de buurt een parkeerplaats vinden. In figuur 10 is te zien hoe het aantal auto’s op het parkeerterrein kan fluctueren in de tijd. Opgave 6. Met de applet Naam kun je de situatie op het parkeerterrein van arriverende en vertrekkende auto’s simuleren, en aan de hand van hiervan schattingen geven voor de verliespercentages. (i) Bepaal de waarden van λ en µ. (ii) Neem aan dat we mikken op een verliespercentage van 0.1%; probeer met behulp van de applet Naam een indruk te krijgen van de gewenste grootte van het parkeerterrein.
2.2
Formule voor het verliespercentage
Waarschijnlijk is het je in de vorige opdracht gelukt om met gewoon proberen een verstandige keus te maken voor de grootte van het parkeerterrein. Maar het is efficienter om de 16
···
−
···
−
Figuur 11: Stroomdiagram van het parkeerterrein met N plaatsen gewenste grootte van het parkeerterrein te bepalen op basis van een wiskundige aanpak. Hieronder zullen we laten zien hoe voor een gegeven parkeerterrein het verliespercentage exact kan worden berekend. De grootte van het parkeerterrein geven we aan met N en de verlieskans (of verliesfractie) met BN ; het verliespercentage is dan BN × 100%. Verder is pn de kans of fractie van de tijd dat er n auto’s op het parkeerterrein zijn (n = 0, 1, . . . , N ). We zijn in het bijzonder geinteresseerd in pN , want λ · pN is het gemiddeld aantal auto’s per uur dat een vol parkeerterrein aantreft, en dus aantal auto’s per uur dat het parkeerterrein vol aantreft totaal aantal arriverende auto’s per uur λ · pN = λ = pN .
BN =
De kansen pn kunnen op intuitieve wijze met behulp van een zogenaamd stroomdiagram worden bepaald. Een stroomdiagram bestaat uit knopen die verbonden zijn met pijlen. De knopen geven de toestanden aan en de pijlen de overgangen die tussen deze toestanden mogelijk zijn. Elke pijl heeft als label de intensiteit van die overgang. Het stroomdiagram voor het parkeerterrein is te zien in figuur 11. De toestanden (aantal auto’s) zijn genummerd van 0 tot en met N . Er zijn twee soorten overgangen mogelijk: aankomsten en vertrekken. Een aankomst van een auto vindt plaats met intensiteit λ en leidt tot een overgang van toestand n naar n + 1. De intensiteit waarmee een auto vertrekt is µ, en dus is de totale vertrekintensiteit in toestand n gelijk aan n · µ (want er staan n auto’s). Een vertrek leidt tot een overgang van n naar n − 1. Voor de kansen pn kunnen we een stelsel vergelijkingen opstellen met behulp van een balansargument: als we het stroomdiagram doormidden snijden tussen toestand n − 1 en n, dan zal gelden dat het aantal sprongen per tijdseenheid van links naar rechts gelijk is aan het aantal sprongen per tijdseenheid van rechts naar links. In figuur 11 zien we dat het alleen vanuit toestand n − 1 mogelijk is om van de linker naar de rechterhelft te springen, en dus geldt het aantal sprongen per tijdseenheid van links naar rechts = λ · pn−1 17
en net zo, het aantal sprongen per tijdseenheid van rechts naar links = nµ · pn . Gelijkstellen van beide stromen levert λ · pn−1 = nµ · pn , ofwel
1 λ · · pn−1 , n = 1, 2, . . . , N. n µ Door herhaald toepassen van deze gelijkheid kunnen we pn als volgt uitdrukken in p0 . 1 λ pn = · · pn−1 n µ 1 λ2 = · 2 · pn−2 n(n − 1) µ = ··· 1 λn = · n · p0 n(n − 1)(n − 2) · · · 2 · 1 µ 1 n = ρ p0 , (9) n! waarin ρ = λ/µ. Verder moet gelden dat de kansen pn optellen tot ´e´en, dus pn =
p0 + p1 + · · · + pN = 1. Als we formule (9) invullen voor pn , krijgen we de volgende vergelijking voor p0 : 1 2 1 N p0 1 + ρ + ρ + · · · + ρ = 1, 2! N! zodat 1 . p0 = 1 2 1 + ρ + 2! ρ + · · · + N1 ! ρN Voor de verliesfractie BN vinden we tenslotte 1 N ρ N! BN = pN = 1 + ρ + 2!1 ρ2 + · · · + N1 ! ρN en dit is precies de beroemde formule van Erlang! Opgave 7. Met behulp van de kansen pn kan de gemiddelde bezetting L van het parkeerterrein worden berekend volgens L = 0 · p0 + 1 · p1 + 2 · p2 + · · · + N · pN . Maar het is eenvoudiger om L te berekenen met de formule van Little. Laat met behulp van de formule van Little zien dat geldt L = ρ · (1 − BN ).
18
2.3
Berekening van het verliespercentage
De formule van Erlang mag er mooi uitzien, maar de formule lijkt niet handig voor numerieke berekeningen. De formule bevat N !. Probeer dat voor grote N (bijvoorbeeld N = 100) maar eens uit te rekenen op een zakrekenmachine! Hieronder zullen we een eenvoudig algoritme af leiden voor de berekening van het verliesfractie BN . Als we in de formule voor BN teller en noemer delen door 1 + ρ + 2!1 ρ2 + · · · + krijgen we 1 N 1 ρ /(1 + ρ + 2!1 ρ2 + · · · + (N −1)! ρN −1 ) N! . BN = 1 ρN −1 ) 1 + N1 ! ρN /(1 + ρ + 2!1 ρ2 + · · · + (N −1)!
1 ρN −1 (N −1)!
Merk op dat 1 N ρ N!
1 + ρ + 2!1 ρ2 + · · · +
1 ρN −1 (N −1)!
1 ρN −1 ρ (N −1)! · 1 N 1 + ρ + 2!1 ρ2 + · · · + (N −1)! ρN −1 ρ = · BN −1 , N
=
en dus geldt voor BN dat BN =
1 N
· ρBN −1 ρBN −1 . = 1 N + ρBN −1 1 + N · ρBN −1
Met deze formule kunnen we BN berekenen zodra we BN −1 hebben berekend, en voor BN −1 moeten we eerst BN −2 berekenen, enzovoort. Met andere woorden, we kunnen BN recursief berekenen als we starten met B0 = 1. Het algoritme is hieronder beschreven. Algoritme voor de berekening van de verliesfractie BN 1. n := 0 B := 1 2. while n < N do n := n + 1 B := rho * B / (n + rho * B) 3. return B De applet Naam gebruikt bovenstaand algoritme voor de berekening van de verliesfractie. Opgave 8.
19
(i) Bepaal met behulp van de applet Naam de gewenste grootte van het parkeerterrein als we mikken op een verliespercentage van 0.1%; vergelijk het antwoord met de resultaten uit Opdracht 6. (ii) Bereken de gemiddelde bezetting van het parkeerterrein gevonden in (i). (iii) Onderzoek of de gewenste grootte van het parkeerterrein gevoelig is voor het verliespercentage, de aankomstintensiteit en de gemiddelde winkeltijd.
Referenties [1] J.A. Buzacott, J.G. Shanthikumar (1993) Stochastic models of manufacturing systems. Englewood Cliffs, Prentice Hall. [2] A.K. Erlang (1909) The Theory of Probabilities and Telephone Conversations. Nyt Tidsskrift for Matematik B, vol 20. [3] A.K. Erlang (1917) Solution of some Problems in the Theory of Probabilities of Significance in Automatic Telephone Exchanges. Elektrotkeknikeren, vol 13. [4] W.J. Hopp, M.L. Spearman (2001) Factory physics: foundations of manufacturing managament. London, McGraw-Hill. [5] N. Gans, G. Koole, A. Mandelbaum (2002) Telephone call centers: a tutorial and literature review. [6] S.S. Lavenberg (1983) Computer performance modeling handbook. London, Academice Press. [7] J.D. Little (1961) A proof of the queueing formula L = λW . Opns. Res., vol 91. [8] H.C. Tijms (1999) Spelen met kansen. Utrecht, Epsilon. [9] H.C. Tijms, F. Heiermans, R. Nobel (2000) Poisson, de Pruisen en de lotto: de Poisson verdeling en haar toepassingen. Utrecht, Epsilon. [10] H.C. Tijms (2002) Operationele analyse: een inleiding in modellen en methoden. Utrecht, Epsilon. [11] J. Walrand, P. Varaiya (2000) High-performance communication networks. San Francisco, Morgan Kaufmann. [12] F.V. Webster (1958) Traffic signal settings. Department of Scientific and Industrial research, Road Research Laboratory, Technical Paper No. 39.
20
