av Øystein Myrland, Professor v/UiT

Vi skal i dette notatet modellere tidsplanlegging under usikkerhet ved å benytte PERT metodologien. Dette gjøres med en Monte Carlo simulering i R. Nødvendige R pakker lastes først.

require(pacman) || {install.packages("pacman"); require(pacman)}
p_load(mc2d, mosaic, reshape2, ggpubr)

I pakken mc2d finner vi funksjoner som genererer PERT fordelingen. Vi skal først se på sannsynlighetsfordelingen (pdf). Selve tetthetsfunksjonen heter dpert(). Her er det tre verdier som må spesifiseres, der min= er optimistisk tid, mode= er mest sannsynlig tid, og max= er pessimistisk tid. I tillegg er det mulig å endre kurtosis til fordelingen, et mål på hvor “tykk” halen i fordelingen er. Dette gjøres ved å endre shape=. I den opprinnelige PERT fordelingen er kurtosis lik 4. Dette er også verdien som benyttes om du ikke spesifiserer en shape verdi. Mer om PERT funksjonene finner du her.

Vi starter med å lage en figur av PERT fordelingen med samme min, mode og max verdi, henholdsvis 4, 5 og 12. Deretter varierer kurtosis mellom 6, 4 og 3. Vi ser at jo mindre kurtosis er, jo “tykkere” blir halen i fordelingen. Dette tolkes som at mer pessimistiske verdier av tidsbruk blir mer sannsynlig.

curve(dpert(x, min=4, mode=5, max=12, shape=6), from = 3, to = 13, main="PERT fordeling med varierende kurtosis", xlab="Tid", ylab="Sannsynlighet", lwd=2)
curve(dpert(x, min=4, mode=5, max=12, shape=4), from = 3, to = 13, add=T, col="blue", lwd=2)
curve(dpert(x, min=4, mode=5, max=12, shape=3), from = 3, to = 13, add=T, col="red", lwd=2)
legend("topright", inset=.05, title="Kurtosis", c("6","4","3"), fill=c("black","blue","red"), horiz=TRUE)

Vi skal nå gjøre en Monte Carlo simulering. Dette gjøres ved å la datamaskinen trekke n tilfeldinge verdier under en spesifikk PERT fordeling. Til dette benyttes funksjonen rpert(). Vi setter en seed verdi som gjør at vi får samme tilfeldige resultat hver gang. Jo flere trekninger en gjør, jo mer “nøyaktig” blir resultatene. En stor grad av nøyaktighet får vi ved å trekker hver aktivitet 10 000 ganger. Dataene henter vi fra Tabell 9.3 Eksempel B, s. 318. Vi velger å benytte oss av en “vanlig” kurtosis på 4. Vi gir aktivitetene navn fra A-J.

set.seed(42)
n=10000 
kurt=4
A <- rpert(n, min=4, mode=5, max=12, shape=kurt)
B <- rpert(n, min=1, mode=1.5, max=5, shape=kurt)
C <- rpert(n, min=2, mode=3, max=4, shape=kurt)
D <- rpert(n, min=3, mode=4, max=11, shape=kurt)
E <- rpert(n, min=2, mode=3, max=4, shape=kurt)
F <- rpert(n, min=1.5, mode=2, max=2.5, shape=kurt)
G <- rpert(n, min=1.5, mode=3, max=4.5, shape=kurt)
H <- rpert(n, min=2.5, mode=3.5, max=7.5, shape=kurt)
I <- rpert(n, min=1.5, mode=2, max=2.5, shape=kurt)
J <- rpert(n, min=1, mode=2, max=3, shape=kurt)

La oss først se på litt deskriptiv statistikk til aktivitet A. Vi ser at minimumsverdien er 4, mens maksimumsverdien er 11.41 for de 10 000 trekningene. Videre så ser vi at gjennomsnittet er 5.99, mens variansen er 1.73.

favstats(A)
      min       Q1   median       Q3      max     mean       sd     n missing
 4.001207 4.940727 5.755889 6.788759 11.40682 5.989334 1.316917 10000       0

En enkelt figur av den simulerte fordelingen til aktivitet A (10 000 observasjoner) kan gjøres med funksjonen densityplot(). Vi benytter denne som et alternativ til et histogram. Jo flere trekninger (n) en gjør i en Monte Carlo simulering, jo “glattere” blir figuren av trekningene. Det er svært enkelt å gjøre en million trekninger, det tar bare litt lenger tid.

densityplot(A, main="Empirisk tetthetsfordeling til aktivitet A", xlab="Tid", ylab="Sannsynlighet", lwd=2)

Nå har vi 10 000 trekninger for hver aktivitet A-J under hver sin PERT fordeling. Deretter genererer vi de fire alternative stiene som er skissert i Figur 9.10 AON-nettverk - eksempel B, s. 319.

ACFJ = A+C+F+J
ADGJ = A+D+G+J
AEHIJ = A+E+H+I+J
BHIJ = B+H+I+J

For å kunne lage noen fine figurer, legger vi de fire stiene i et datasett, kalt “stier”.

stier = melt(cbind(ACFJ,ADGJ,AEHIJ,BHIJ))
names(stier) = c("N","Sti","Tid")
head(stier)
  N  Sti      Tid
1 1 ACFJ 12.32798
2 2 ACFJ 13.25044
3 3 ACFJ 13.39705
4 4 ACFJ 14.37998
5 5 ACFJ 12.93124
6 6 ACFJ 17.35275

Vi lager deretter en figur av de fire stiene sin sannsynlighetstetthet (pdf).

ggdensity(stier, x = "Tid", add = "mean", color = "Sti", fill = "Sti") + labs(title = "Sannsynlighetstettheter (pdf)") + ylab("Sannsynlighet")

Vi kan også lage en figur av den kumulative sannsynlighetstettheten (cdf). I denne figuren er det enkelt å se hvilken sti som er den kritiske. Den kritiske stien vil alltid ligge lengst til høyre i en slik fremstilling. Den har alltid høyest timetall med en gitt sannsynlighet. Sammenligner en svært mange stier er et cdf plott den enkleste metoden for å finne kritisk sti.

ggplot(stier, aes(x=Tid, colour = Sti)) + stat_ecdf() + labs(title = "Kumulativ Sannsynlighet (cdf)") + ylab("Sannsynlighet")

Den kumulative sannsyligheten gir oss sannsynligheten for at \(P(Tid \le x)\). Vi velger et timetall (Tid) på x-aksen, går deretter opp til vi treffer den kumulative sannsynligheten til en Sti, og leser av sannsynligheten av på venstre akse som \(P(Tid \le x)\).

Det blir lett unøyaktig dersom vi skal lese dette av grafisk. Vi kan derfor generere en empirisk kumulativ sannsynlighetsfunksjon for de fire stiene, slik at vi kan regne på sannsynlighetene.

stiACFJ=ecdf(ACFJ)
stiADGJ=ecdf(ADGJ)
stiAEHIJ=ecdf(AEHIJ)
stiBHIJ=ecdf(BHIJ)

Når en har disse kan en beregne spesifikke sannsynligheter for hver sti. I eksempelet er kritisk sti AEHIJ. Sannsynligheten for at denne aktiviteten kan fullføres innen 20 tidsenheter, dvs. \(P(X \le 20)\) er:

stiAEHIJ(20)
[1] 0.95

Det betyr at det er 95% sannsynlighet for at aktiviteten kan gjennomføres på 20 eller færre tidsenheter. I boka på s. 319 er denne sannsynligheten beregnet til 96.6%, men vår beregning er mer nøyaktig. Hvordan kan jeg si det? I boka standardiserer en verdien 20 med å benytte standardavvik og gjennomsnitt. En går da fra en positiv skjevfordelt PERT fordeling med en kurtosis på 4, til en symmetrisk standard normalfordeling. Dette blir litt unøyaktig sammenlignet med vår metode som gjør beregninger på selve PERT fordelingen.

Vi kan også se på sannsynligheten for at \(P(X \gt 18)\), denne er \(1-P(X \le 18)\), som er:

1-stiAEHIJ(18)
[1] 0.2554

eller 25.54%. Vi kan også finne sannsynligheten for at \(P(16 \le X \le 20)\), med:

stiAEHIJ(20)-stiAEHIJ(16)
[1] 0.6472

Vi avslutter med å se på litt deskriptiv statistikk til de ulike stiene.

favstats(~Tid, groups = Sti, data = stier)
    Sti       min        Q1   median       Q3      max      mean       sd     n missing
1  ACFJ  9.725866 11.944300 12.80460 13.88916 19.13269 12.992706 1.434582 10000       0
2  ADGJ 11.058508 14.549364 15.78678 17.24421 25.58422 15.972837 1.974009 10000       0
3 AEHIJ 12.622855 15.762066 16.83033 18.03419 23.72921 16.985780 1.676151 10000       0
4  BHIJ  6.977803  9.146123  9.90666 10.77049 14.62731  9.997234 1.165188 10000       0

Gitt denne R koden, så kan du selv raskt se hvor følsom den kritiske stien er for endringer i kurtosis fra 4 til 3, dvs. å gjøre mer pessimistiske verdier for aktivitetene mer sannsynlig. Det er også enkelt å endre koden slik at den passer til et nytt case med fler eller færre aktiviteter, og varierende min, mode, max og shape verdier.

R koden uten kommentarer er som følger:

require(pacman) || {install.packages("pacman"); require(pacman)}
p_load(mc2d, mosaic, reshape2, ggpubr)
curve(dpert(x, min=4, mode=5, max=12, shape=6), from = 3, to = 13, main="PERT fordeling med varierende kurtosis", xlab="Tid", ylab="Sannsynlighet", lwd=2)
curve(dpert(x, min=4, mode=5, max=12, shape=4), from = 3, to = 13, add=T, col="blue", lwd=2)
curve(dpert(x, min=4, mode=5, max=12, shape=3), from = 3, to = 13, add=T, col="red", lwd=2)
legend("topright", inset=.05, title="Kurtosis", c("6","4","3"), fill=c("black","blue","red"), horiz=TRUE)
set.seed(42)
n=10000 
kurt=4
A <- rpert(n, min=4, mode=5, max=12, shape=kurt)
B <- rpert(n, min=1, mode=1.5, max=5, shape=kurt)
C <- rpert(n, min=2, mode=3, max=4, shape=kurt)
D <- rpert(n, min=3, mode=4, max=11, shape=kurt)
E <- rpert(n, min=2, mode=3, max=4, shape=kurt)
F <- rpert(n, min=1.5, mode=2, max=2.5, shape=kurt)
G <- rpert(n, min=1.5, mode=3, max=4.5, shape=kurt)
H <- rpert(n, min=2.5, mode=3.5, max=7.5, shape=kurt)
I <- rpert(n, min=1.5, mode=2, max=2.5, shape=kurt)
J <- rpert(n, min=1, mode=2, max=3, shape=kurt)
favstats(A)
densityplot(A, main="Empirisk tetthetsfordeling til aktivitet A", xlab="Tid", ylab="Sannsynlighet", lwd=2)
ACFJ = A+C+F+J
ADGJ = A+D+G+J
AEHIJ = A+E+H+I+J
BHIJ = B+H+I+J
stier = melt(cbind(ACFJ,ADGJ,AEHIJ,BHIJ))
names(stier) = c("N","Sti","Tid")
head(stier)
ggdensity(stier, x = "Tid", add = "mean", color = "Sti", fill = "Sti") + labs(title = "Sannsynlighetstettheter (pdf)") + ylab("Sannsynlighet")
ggplot(stier, aes(x=Tid, colour = Sti)) + stat_ecdf() + labs(title = "Kumulativ Sannsynlighet (cdf)") + ylab("Sannsynlighet")
stiACFJ=ecdf(ACFJ)
stiADGJ=ecdf(ADGJ)
stiAEHIJ=ecdf(AEHIJ)
stiBHIJ=ecdf(BHIJ)
stiAEHIJ(20)
1-stiAEHIJ(18)
stiAEHIJ(20)-stiAEHIJ(16)
favstats(~Tid, groups = Sti, data = stier)

Denne R koden kan du klippe ut, lime inn i http://rstudio.uit.no, og kjøre den der. Logg deg på med hele epostadressen din, inklusiv alfakrøll.uit.no. Er du utenfor campus, husk å starte VPN først.