Probabilistic critical path analysis - a comprehensive survey

PROBABILISTIC CRITICAL PATH ANALYSIS A COMPREHENSIVE SURVEY Stefan Isermann

CONTENTS 1. Critical Path Analysis – Deterministic Methods 2. Uncertainty and probabilities and how they affect project duration 3. Main results from statistical analysis 4. Criticality and severity of activities 5. Algorithms A Mathematical Appendix References

Stefan Isermann/ Critical Path Analysis

27.08.2015

2

1. Critical Path Analysis – Deterministic Methods Introduction Probably most people dealing with projects have heard of the critical path in a project and I suppose many of them have got an idea of the meaning of the critical path. I am convinced that all who are certified as a Project Management Professional (PMP) are able to find the deterministic critical path and to calculate project duration in a given project plan. But I am convinced too that are only a few project managers who really understand the probabilistic interaction of paths and activities which leads to a longer project duration than calculated by your project management tool. To deliver just in time (and quality) is one of the main challenges in complex project management. So every project manager should be aware of the risk of deterministic project planning. The aim of this presentation is to enable project manager to understand which calculations or simulations are necessary to get probabilistic statements about project duration and furthermore make proposals how to implement appropriate algorithms.

Stefan Isermann / Critical Path Analysis

27.08.2015

3

1. Critical Path Analysis – Deterministic Methods Definitions In a project management course you may have learned how to calculate the critical path. Let us take a short review.

The critical path method calculates the minimum completion time for a project. The critical path itself represents the set or sequence of predecessor/successor activities which will take the longest time to complete. Forward pass: Earliest event times are computed as the maximum of the earliest start times plus activity durations for each of the activities immediately preceding an event. The earliest start time for each activity (i,j) is equal to the earliest possible time for the preceding event E(i) Backward pass: The procedure for finding the latest event time is analogous to that for the earliest event time except that the procedure begins with the final event and works backwards through the project activities. Critical path: Events which have equal earliest and latest start (or end) time lie on the critical path. Let us take a look at the following example to get a clear understanding. It is an AON (activity on node) representation of a simpler looking AOA (activity on arc) model in pmbook (chapter 10, Fig. 10-4)

Stefan Isermann / Critical Path Analysis

27.08.2015

4

1. Critical Path Analysis – Deterministic Methods Find critical path by forward and backward passing ES, EF denote the earliest start and finish time and LS, LF the latest. Project duration and activity name is placed in the middle. Let us do the forward pass And now the backward pass Activities on the critical path are those with EF=LF and ES=LS. The calculated project duration is 30.

Stefan Isermann / Critical Path Analysis

27.08.2015

5

1. Critical Path Analysis – Deterministic Methods Another heuristic method A more heuristic method is checking all passes and identify the longest as the critical path. This method seems not to be as elegant as the previous algorithm. And do we have really checked all possible paths?

Stefan Isermann / Critical Path Analysis

27.08.2015

6

1. Critical Path Analysis – Deterministic Methods Advantage of the heuristic method Although the method seems to be laborious it provides some insight into the underlying project structure. You see that there is another (green) path whose duration is very near to the critical path duration and obviously may affects or becomes the critical path if there is some delay in activity E or G. If the duration of activities becomes probabilistic on the way from plan to reality, there might be an interaction of all possible path. In this example we have identified 8 different paths just by trying out. To study a deterministic method which contributes the number of paths as well as the critical path length we take a short trip to graph theory.

Stefan Isermann / Critical Path Analysis

27.08.2015

7

1. Critical Path Analysis – Deterministic Methods Topological Sorting of a Directed Acyclic Graph Each project plan may be represented by a Directed Acyclic Graph (DAG). That means, it is formed by a collection of nodes and arrows (directed), where there is no way to start from one node following a sequence of arrows and come back to the node again (acyclic). A topological ordering of a directed graph is a linear order where each node u which is followed by a node v comes before v in the ordering. Here you can see a possible topological order of our project. If you are able to find a topological sorting of a graph it has to be a DAG. And vice versa if you have a DAG you can always do a topological sorting, whereas the topological order is not unique. See another example. Stefan Isermann / Critical Path Analysis

27.08.2015

8

1. Critical Path Analysis – Deterministic Methods Critical path length and number of paths of a topological sorted DAG Once you have a topological sorted representation of your project the length of the critical path is straight forward. At each node u the length of the critical path lu is the duration of the node du plus the maximum duration of all predecessors: lu= du + max(lpred(u)) So you get a critical path length of 30. And by simple backtracking, taking the maximum of the predecessor into account, the path itself: I-F-C-A Whereas the number of paths at each node nu is just the sum of paths of all predecessors : nu= Σ npred(u) So the number of 8 paths is confirmed.

node

d

S

0

A

4

B

pred

l

n

0

1

S

4+0=0

1

3

S

3+0=3

1

C

8

A

8+4=12

1

D

7

A

7+4=11

1

E

9

B, C

9+12=21

1+1=2

F

12

B, C

12+12=24

1+1=2

G

2

D, E

2+21=23

1+2=3

H

5

D, E

5+21=26

1+2=3

I

6

F, G

6+24=30

2+3=5

End

0

H, I

0+30=30

3+5=8

Stefan Isermann / Critical Path Analysis

27.08.2015

9

1. Critical Path Analysis – Deterministic Methods Some remarks on topological sorted DAG Topological sorting algorithm for DAG are very simple and easy to implement. Furthermore they are very quick because of linear order especially O(|nodes| + |edges|). A topological sorted DAG may serve as an input for some critical path algorithms especially for simulations as for algorithms we discuss below. Unfortunately it is not possible to study the interactions of all paths simultaneously. In the previous example the number of path is manageable but in large projects of normal complexity, the number of paths increases exponentially and may exceed the number of atoms in our observable universe! If you do not believe see App. A.8.

100.000.000.000.000.000.000. 000.000.000.000.000.000.000. 000.000.000.000.000.000.000. 000.000.000.000.000.000 The detailed, all-sky picture of the infant universe created from nine years of WMAP data. The image reveals 13.77 billion year old temperature fluctuations (shown as color differences) that correspond to the seeds that grew to become the galaxies. This image shows a temperature range of ± 200 microKelvin. (Credit: NASA / WMAP Science Team)

Stefan Isermann / Critical Path Analysis

27.08.2015

10

CONTENTS 1. Critical Path Analysis – Deterministic Methods 2. Uncertainty and probabilities and how they affect project duration 3. Main results from statistical analysis 4. Criticality and severity of activities 5. Algorithms A Mathematical Appendix References

Stefan Isermann/ Critical Path Analysis

27.08.2015

11

2. Uncertainty and how they affect project duration Uncertainty leads to probabilities If we accept that there is some uncertainty in estimation of activity duration, we have to ask what is the probability of another path (the green one) to become the critical one and what happens with the expected project duration? At first you can treat path duration as a random variable obeying any probability distribution. According to the distributions sometimes the green path is longer than the red path and vice versa (at this moment we neglect the other paths). As the critical path randomly changes according two the maximum of the two paths, we like to know the probability distribution of the maximum and particularly the expected value of the maximum which equals the expected project duration. Since not everyone is familiar with statistical methods, we like to illustrate the problem by a small example. Let us throw two dice and ask the question: “What is the probability that the maximum of the two dice is less or equal 4”? Do you know the answer? First think and then go to the next slide…

Stefan Isermann / Critical Path Analysis

27.08.2015

12

2. Uncertainty and how they affect project duration Throw dice max

1 2 3 4 5 6

1 1 2 3 4 5 6

2 2 2 3 4 5 6

3 3 3 3 4 5 6

4 4 4 4 4 5 6

5 5 5 5 5 5 6

6 6 6 6 6 6 6

cdf of maximum of two dice 100% 80% 60% 40% 20% 0%

You may just count positive events and divide it by all events: P(max(d1,d2) ≤ 4) = 16/36 = 4/9 But obviously the probability that the maximum of the two dice is smaller or equal 4 is the probability that dice 1 and dice 2 is smaller or equal 4! P(max(d1,d2) ≤ 4) = P(d1 ≤ 4) P(d2 ≤ 4) = 4/6 ∙ 4/6 = 4/9 P(max(d1,d2) ≤ n) = n2/36 is called the cumulative distribution function (cdf) of the maximum of two dice Theorem: Let X1, …, Xn be independent random variables each with cdf Fi(x) = Pi(Xi ≤ x) and Xmax = max(X1, …, Xn ) then P(X … ∙FFn(x) (2.1) P(Xmax ≤ x) = FX (x) = F1(x) ∙FF2(x) ∙… The cumulative distribution function of the maximum on n independent random variables is the product of the cumulative distribution function of each variable! max

1

2

3

4

Fig. 2.1

5

6

Stefan Isermann / Critical Path Analysis

27.08.2015

13

2. Uncertainty and how they affect project duration Normal distribution Now you have got the key point to all further considerations. Let us assume, that path duration is continuously and normally distributed with mean µ and standard deviationσ and denote the standard normal probability density function (pdf (pdf) pdf) ϕ(x) and the corresponding cumulative distribution function (cdf (cdf) cdf) Φ(x) For a continuous distribution the pdf (probability density function) f (x) is defined as the derivative of the cdf dF ( x) (2.2) f ( x) =

1 − x2 / 2 e , 2π 1  x−µ  f ( x, µ , σ ) = ϕ   σ  σ 

ϕ ( x) =

dx

x

Φ ( x) = ∫ ϕ ( z ) dz = −∞

1 2π

 x−µ  F ( x, µ , σ ) = Φ    σ 

x

∫e

−∞

−z2 / 2

dz

If we assume that paths are independent and duration is normally distributed the calculation for the cdf and pdf of project duration is straight forward.

Fig. 2.2 Stefan Isermann / Critical Path Analysis

27.08.2015

14

2. Uncertainty and how they affect project duration Probability distribution of maximum of independent normal variables With (2.1) P(Xmax ≤ x) = FX (x) = F1(x) ·F2(x) ·… ·Fn(x) we get the cdf of the maximum of n independent normal distributed random variables  x−µ  (2.3)  F ( x ) = ∏Φ  max

n

i

i =1

 σi 

To simplify matters let us assume that we have only two independent random variables. Then we get for cdf F(x) and pdf f (x) of the maximum  x−µ   x−µ  Φ   F ( x ) = Φ  (2.4)  σ   σ  1

1

f ( x) =

2

2

dF ( x) 1  x − µ1   x − µ 2  1  x − µ1   x − µ 2  Φ  + Φ ϕ   = ϕ  dx σ 1  σ 1   σ 2  σ 2  σ 1   σ 2 

(2.5)

To better understand eq. (2.5) we show the pdf f (x) of the maximum of two independent normal variables for different parameters σi, µ i. Stefan Isermann / Critical Path Analysis

27.08.2015

15

2. Uncertainty and how they affect project duration Typical distributions of maximum

Fig. 2.3a

Fig. 2.3b

Fig. 2.3c

Fig. 2.3d

It is obvious that the distribution of the maximum depends on the ‘distance’ of the respective distributions. In the first and last figure, the distribution of the maximum (the project duration) is nearly the distribution of the variable with greatest mean. If the distributions are adjacent (fig. 2.3b and 2.3c), the distribution of the maximum is right skewed and does not show a Gaussian shape and the mean is shifted to larger values. In the next chapter we will quantify the meaning of distance, adjacent and shift. Stefan Isermann / Critical Path Analysis

27.08.2015

16

2. Uncertainty and how they affect project duration Assumption of normal distribution and independence Up to this point one may suppose that the problem of critical path analysis is solved by. eq. (2.1) in principle. But there are two questions left: 1. Is the duration of a path normally distributed if the duration of an activity is estimated via PERT (beta distribution with shape parameters α +β = 6)? 2. Is it correct to assume durations of paths as independent (not correlated) random variables? The answers: 1. Due to central limit theorem, the fact that the shape of a symmetric beta distribution (PERT) is close to a normal distribution and the maximal skewness γ is limited to|γ|< (7/5)½ = 1.18 the distribution of path duration is close to a normal distribution even for a small number of activities! 2. Unfortunately and in general different paths in a project plan have some activities in common. E. g. the red and the green path in our sample project have activity A, C, G and I in common. So paths in a project are generally dependent and correlated and the simple eq. (2.3) can not be applied.

Stefan Isermann / Critical Path Analysis

27.08.2015

17

2. Uncertainty and how they affect project duration Normal distribution is confirmed, independence not To confirm the first statement let us take a look at the path distributions of our sample project, assuming that all activities are beta distributed. As an example all activities have a relative standard deviation of 25% with a right skewed beta distribution (α = 2, β = 4, γ = 0.47) The figures (obtained via convolution integral or fast fourier transformation, s. App. A.1) shows that even for the small number of activities in our sample project the path durations exhibit nearly normal distributions. So we are left with the problem of path correlation.

Fig. 2.4 Stefan Isermann / Critical Path Analysis

27.08.2015

18

CONTENTS 1. Critical Path Analysis – Deterministic Methods 2. Uncertainty and probabilities and how they affect project duration 3. Main results from statistical analysis 4. Criticality and severity of activities 5. Algorithms A Mathematical Appendix References

Stefan Isermann/ Critical Path Analysis

27.08.2015

19

3. Main results from statistical analysis Exact cumulative distribution function of maximum Each path in a project may have correlations to other paths in a project depending on the number of common activities and their variances. It is straight forward to calculate the corresponding covariance matrix Σ of all project paths (s. App. A.2). As shown before we may assume that the common distribution of all k paths in a project is the multivariate normal distribution with density 1 (3.1)  1  f (x , x ,L, x ) = exp − (x − µ ) Σ (x − µ ) T

1

2

k

(2π) Σ

−1

 2

k



The probability (cdf) that the maximum of all paths is smaller than z is z

z

z

−∞

−∞

−∞

F(max(x1 , x 2 , L , x k ) = z ) = ∫ dx1 ∫ dx 2 L ∫ dx k

 1  T exp − (x − µ ) Σ −1 (x − µ ) k 2   ( 2 π) Σ 1

(3.2)

If Σ is diagonal (no correlations) the integrand in (3.2) factorizes and we directly get eq. (2.3). Unfortunately in all other cases the integral (3.2) can be solved only numerically even in two dimensions! Stefan Isermann / Critical Path Analysis

27.08.2015

20

3. Main results from statistical analysis Exact results in two dimensions But in two dimensions the pdf of the maximum of two correlated Gaussian variables may be expressed in terms of the normal probability and cumulative function. All moments of the maximum follows by integration. You will find all these results in App. A.3. These results are first obtained by C. E. Clark [1]. Some of these results are provided and discussed in a short paper of Nadarajah, Kotz [7]. The expectation value of the maximum of two Gaussian random variables z = max(x1, x2) is µ −µ  µ −µ  µ −µ  E( z) = µ Φ  + µ Φ  + θϕ   (3.3)  θ   θ   θ  with θ = σ + σ − 2 ρσ σ (3.4) If we define a normalized distance α between the distributions of the two Gaussian random variables µ −µ α= (3.5) θ we can rewrite eq. (3.3) E (z ) = µ + θ (ϕ (α ) − αΦ (− α )) (3.6) 1

2

2

1

2 1

2 2

1

1

1

2

2

1

2

2

1

Stefan Isermann / Critical Path Analysis

27.08.2015

21

3. Main results from statistical analysis Shift of maximum of two Gaussian random variables Without loss of generality we assume µ1 ≥ µ2 (α ≥ 0). So the addend in eq. (3.6) is just the shift of the maximum. The shift of the maximum to higher values is larger if correlation between paths decreases (θ increases) and the normalized distance between paths decreases (α decreases) and is largest if µ2 = µ1 (α =0). The maximal shift of the maximum of two Gaussian random variables is θ shift = ≅ 0.4 ⋅ θ (3.7) 2π And with µ2 = µ1 and σ1= σ2= σrel µ1 we get the order of the maximal relative shift σ rel. shift = ≅ 0.56 ⋅ σ (3.8) π So depending on the amount of uncertainty and the project structure the shift of project duration may be of the order of some percent up to 10%. Fig. 3.1 max

rel

max

rel

Stefan Isermann / Critical Path Analysis

27.08.2015

22

3. Main results from statistical analysis Comparison of the shift of the maximum Now we come back to the figures of chapter 2. Fig. 3.2a and 3.2d correspond to a normalized distance between two Gaussians with α = 2. The shift is small and not visible to the naked eye. Furthermore the pdf of the maximum is almost Gaussian. Fig. 3.2b and 3.2c correspond to a normalized distance α = 0.5. The shift is visible and about 20 times larger. The distribution in fig. 3.2b and 3.2c are obviously right skewed, whereas in 3.2b it is more right skewed than in 3.2c. This is due to the fact that skewness is not symmetric in α and larger for negative α if the index 1 in α is coupled to the Gaussian with greater standard deviation σ (s. App. A.4).

Fig. 3.2a

Fig. 3.2b

Fig. 3.2c

Fig. 3.2d

Stefan Isermann / Critical Path Analysis

27.08.2015

23

CONTENTS 1. Critical Path Analysis – Deterministic Methods 2. Uncertainty and probabilities and how they affect project duration 3. Main results from statistical analysis 4. Criticality and severity of activities 5. Algorithms A Mathematical Appendix References

Stefan Isermann/ Critical Path Analysis

27.08.2015

24

4. Criticality and Sensitivity of Activities Definitions The criticality index (CI) of each activity is defined as the probability that each activity is on a critical path of the project CI =P (CriticalPath ) (4.1) To identify those activities whose uncertainty in duration mostly affects project duration you may calculate the schedule sensitivity index (SSI), which measures the relative importance of an activity sd CI SSI = (4.2) sd where sd means the standard deviation of the respective duration. To get the CI of an activity you have to sum up all probabilities of paths where the activity is member of. For these and all other sensitivity measures you need the probability that a path j is a critical path 1  1  (4.3) P (x > max( x , x ,L , x , x ,L , x ) ) = dx dx dx L dx dx L dx exp − (x − µ ) Σ (x − µ ) activity

activity

activity

activity

activity

project

∞

j

1

2

j −1

j +1

k

∫

−∞

xj

j

∫

−∞

xj

1

∫

−∞

xj

2

∫

−∞

xj

j −1

∫

−∞

xj

j +1

∫

−∞

T

k

(2π ) Σ k

 2

−1



Stefan Isermann / Critical Path Analysis

27.08.2015

25

4. Criticality and Severity of Activities Path probability Once more the integral (4.3) can analytically and easily be solved for two Gaussian random variables by the linear transformation z = x1- x2  1  z − µ1 + µ 2  2  1 P(x1 > x2 ) = ∫ dz exp −    = Φ (α ),  2 θ 2 π θ   0  µ −µ θ = σ 12 + σ 22 − 2 ρσ 1σ 2 , α = 1 2 ∞

(4.4)

θ

For two Gaussian distributed paths the probability to be the critical path is completely determined by the normalized distance α. For |α| > 4 the probability for interaction is smaller than 0.003% so that the only path with greater mean will survive.

Fig. 4.1 Stefan Isermann / Critical Path Analysis

27.08.2015

26

CONTENTS 1. Critical Path Analysis – Deterministic Methods 2. Uncertainty and probabilities and how they affect project duration 3. Main results from statistical analysis 4. Criticality and severity of activities 5. Algorithms A Mathematical Appendix References

Stefan Isermann/ Critical Path Analysis

27.08.2015

27

5. Algorithms Introductory remarks You may find many papers (see references in [9]) concerning the research on probabilistic project completion time by critical path analysis. Unfortunately surfing the web you will find a lot of naive suggestions. Some advise to add just the variances (calculated from PERT) of the deterministic CP. Other argue that independence of activities implies independence of path. All this is wrong. As we have seen the real world looks more complicated. Adding activity duration to different paths will correlate these paths. Repeating this process will end up in a great number of highly correlated paths, which are candidates for a critical path (means high probability to be a critical path). And the number of paths in a project increases exponentially with the number of activities. So it is all but impossible to take into account all or only the most critical paths (except for very small projects in the range of 10 - 20 activities) . Topological sorting is a way out for all promising algorithms. The algorithm for topological sorting is well known, consists of a few lines of code and has linear running time. So we start with a topological sorted DAG (Directed Acylic Graph) as an input for all considered algorithms.

Stefan Isermann / Critical Path Analysis

27.08.2015

28

5. Algorithms Introductory remarks To keep it simple we only consider activity interactions with predecessor finish to successor start constraints (FS). All other activity interactions may be modeled by introducing dummy activities of appropriate duration (including negative values) with zero variance. First example: activity A should begin 5 days after the start of activity B . Second example: activity A should finish 5 days before the end of activity B.

5D2

A

A

-5D4

0D1

B

B

0D3

To write down all rules for lead and lag constraints is not the aim of this paper. Algorithm to implement such rules are straight forward. And at least: every analytical approach should have a shorter running time than a simple Monte Carlo Simulation (MCS) of the same precision.

Stefan Isermann / Critical Path Analysis

27.08.2015

29

5. Algorithms Monte Carlo Simulation – number of simulations Not only easy to implement, but also necessary in order to study general behavior of distributions of project completion time and to compare with other approaches, we take a closer look at MCS. It is well known that the relative error of the standard deviation σs of a sample of size N from a normal distribution is approximately 1/(2N)1/2 *. So if we choose N=10.000 we get high confidence that the relative error of the estimated standard deviation is less than 1%. With σs estimated from the sample the relative error of the estimated mean is σs σµ = which is typically of the order of 0.1%. µ N For a sample resulting from normal distribution the standard error of estimated skewness is approximately (6/N)1/2. So if we estimate skewness from MCS with N=10.000, we may state that the distribution does not differ significantly from normal if skewness is smaller than 0.05. As a consequence we choose N=10.000 in all simulations, because this should be sufficient for all project estimations. All following simulations and numerical calculations are done in R [11]. R is a language and environment for statistical computing and graphics and available as free software. * It differs for all other distribution by an amount of kurtosis, e.g. in case of a symmetric Pert with negative kurtosis γ2 = -2/3 we get sd(σs ) = 1/(3N)1/2

Stefan Isermann / Critical Path Analysis

27.08.2015

30

5. Algorithms Monte Carlo Simulation – simulation pseudo code With just a few lines (10 – 20) of R code you will get a MCS Input: topological sorted adjacency list, vector of means,

vector of standard deviations or covariance matrix for normal dist. or shape and minimum, maximum parameters for beta

Output: mean,

standard deviation, skewness of project duration, probabilities of activities to be on a critical path

Algorithm:

generate ni random values for each activity i with specific distribution parameters from input (to generate correlated skewed random values you need the technique of copulas!) for each activityi from end to start { if number(successorsi) > 1 { activityi,j=1 to N = activityi,j=1 to N + max(activitysuccessorsi,j=1 to N ) compute probability of each successori to be the maximum update probabilities of contributing successor activities else activityi,j=1 to N = activityi,j=1 to N + activitysuccessori,j=1 to N add probability 1 for activity i } }

Stefan Isermann / Critical Path Analysis

27.08.2015

31

5. Algorithms Monte Carlo Simulation – simulated project With just a few lines of code you get this result for our previous example. MC Simulation is done with N=10000 symmetric beta distributed (α=β=3) random variables with mean µi as in chapter 1 and a standard deviation sdi = 0.25µi. The upper figure shows a shift of project completion time from 30 to 31.1, whereas the shape of the distribution is completely normal, which is indicated by the vertical green (mean) and red line (standard deviation). The corresponding p-values of normal distribution (15.9%, 50%, 84.1% - horizontal lines) coincide perfectly with the cumulative distribution of the simulation. The second figure shows that activity E and G which are not part of the deterministic critical path have a great probability (about 40%) to lie on a critical path. Fig. 5.1 Stefan Isermann / Critical Path Analysis

27.08.2015

32

5. Algorithms Monte Carlo Simulation – project simulation with different projects Now we do successive simulations with normal, normal correlated (with ρ = 0.5) and beta right skewed (with constant γ = 0.47, α = 2, β = 4) distributed durations with same mean and standard deviation as above. Whereas the first figure shows no difference to Fig. 5.1 we can see by the second figure that correlation will increase standard deviation and decrease the shift as it is expected and explained in appendix A.2. Furthermore correlation will generally promote activities (increase probability) of the main critical path. The third figure shows that skewness only has a small effect on standard deviation and negligible effect on mean. The deviation from normal is small. Skewness of 0.22 means a deviation of about 1% in probability of mean (γ /(72π)1/2 =1.5%) which we may estimate from the Edgeworth expansion (eq. A4.2). Fig. 5.2 Stefan Isermann / Critical Path Analysis

27.08.2015

33

5. Algorithms Monte Carlo Simulation – correlation and skewness Because expected non-normality in distributions for project completion time is the main point of criticism against analytical approaches, we consider the evolution of skewness with increasing number of activities. Each data point is generated by a simulation of 20 random projects with fixed number of activities as described in Appendix A.7. Complexity of 1.9 means the number of randomly generated edges is approximately 1.9 times the number of activities. The first and third figure show that distribution of completion time of normally distributed (generally symmetric) and even correlated activities in our precision is indistinguishable from normal. And even for projects with activities all right skewed (here γ = 0.47) skewness decreases by approx. 1/N ½ (red line). Fig. 5.3 Stefan Isermann / Critical Path Analysis

27.08.2015

34

5. Algorithms Monte Carlo Simulation - summary To sum up, MC simulations show that the distribution of project completion time of symmetric distributed activities is practically a normal one even for a small number of activities or correlated symmetric activities. Skewness (and therefore all cumulants of higher order) doesn’t matter if number of activities are in the order of 100, because it diminishes as we may expect from central limit theorem in particular the Berry-Esseen theorem approximately by 1/N ½, whereas it is clear that the theorem does only hold for independent random variables. The case of skewed and correlated activities is not a trivial one, as we know that only normally distributed variables are fully described by mean and covariance (i.e. strength of linear correlation), whereas the joined distribution of not normally distributed and dependent variables is described by their marginal distributions and a Copula which describes the structure of dependencies of the variables. That means that you have to know the strength as well as the special kind of dependency structure (a Gaussian copula of defined correlation will lead to different linear correlations coefficients if the marginal distributions are not normal). So we leave this sophisticated discussion to Appendix A.8 whereas it is totally unclear from where you will get information about structure and strength of dependencies of future activities.

Stefan Isermann / Critical Path Analysis

27.08.2015

35

5. Algorithms Analytical approximation In the early sixties, where computational power was far beyond today, Charles E. Clark suggested an algorithm based on the analytical solution of the maximum of two Gaussian random variables [1]. In the meantime many other approaches have been made arguing that the assumptions of normality may be a very poor approximation [7]. Indeed, in the case where two Gaussian variables are strongly neighbored (in the sense of small normalized distance α), and their standard deviation differs by a factor of 3 or greater, the resulting distribution may become strongly skewed. Whereas such cases in project planning are rather pathological we may see in A.4 that for a great region of possible interaction the skewness of the maximum distribution is small (less than 0.5 in skewness means a deviation in probability of the mean of less than 3%). Furthermore it is shown in A.5 by a bivariate Edgeworth expansion that the effect of skewness on the mean of the distribution is negligible and on variance is only a few percent even for strongly skewed beta distributions. To minimize the effect of skewness I adapt an algorithm suggested by Sinha et al. [8], which suggest to minimize the error in maximum operations by a pair-wise max operation which minimizes the integrated absolute difference between exact distribution and normal approximation. Stefan Isermann / Critical Path Analysis

27.08.2015

36

5. Algorithms Analytical approximation - Algorithm It is shown in appendix A.4 and A.5 that the error of maximum approximation is primarily controlled by skewness. As skewness follows directly from eq. A3.4 - 3.6, we try to minimize the approximation error by combining successively those distributions which will have smallest skewness (which is always greater than 0), provided that are more than two distributions (successors). Input of the algorithm is the same as for MCS: adjacency list of topological sorted activities, vector of means, covariance matrix. Go from end to start (if adjacency list is a list of predecessor from start to end), do add – operations for mean and covariance as usual (see appendix A.1), do max operations with the help of formulas A3.1 – 3.6 to compute mean, variance and skewness of the resulting distribution. If you have to compute the maximum of more than two successors, do it pair wise by minimizing the resulting skewness. After each operation update covariance matrix with the help of formulas A6.2 – 6.4. For max – operations update probabilities of involved activities according to A3.3. For each knot you will get the distribution of duration from knot to end as a normal approximation and at last for the whole project. Stefan Isermann / Critical Path Analysis

27.08.2015

37

5. Algorithms Analytical approximation – example project Of course you need some more lines of coding (≈ 80) than for simple MCS but the results are very close to the results of the MCS in fig. 5.2. Even they coincide within the precision of the MCS for 1% in the standard deviation and ± 0.04 or ± 0.06 (ρ = 0.5) for the mean. Furthermore this is the same approximation for skewed distributed activities, which is here close to the normal as we have seen in fig. 5.2. Now we have to answer the question if the goodness of the analytical approach is just a feature of this example project or a general one.

Fig. 5.4 Stefan Isermann / Critical Path Analysis

27.08.2015

38

5. Algorithms Analytical approximation – in general Once more for a different number of activities we generate 20 random projects of complexity cg=1.9 respectively and compare the results of analytical approximation to those of MCS where we compute the mean (points) of the relative error of deviation from MCS - mean and MCS - standard deviation and the corresponding standard deviation (lines). The relative error of approximated mean is of the order of 0.1% regardless of correlation or skewness of all activities. The relative error of approximated standard deviation is less than 1% except for right skewed activities, where the approximation is about 2% less than calculated from MCS. This is what we may expect from our analytical approximation. Overall the result is reasonably good. But what is about computation time? Fig. 5.5 Stefan Isermann / Critical Path Analysis

27.08.2015

39

5. Algorithms Computation time of MCS vs. analytical approximation

Fig. 5.6

Fig. 5.7 Stefan Isermann / Critical Path Analysis

27.08.2015

40

5. Algorithms Computation time of MCS vs. analytical approximation When I do the first MCS of large projects (about 1000 activities) they last about 1 minute till I red the paper of Tattoni et. al. [10] (“Monte Carlo Simulation strikes back”). Encouraged by this paper I improved the code (carefully vectorization in R) up to that what is shown in fig. 5.6. Time measurement is done on a standard NB (Intel Core i5 – 3320M / 2x 2,6 GHz / 2x 4GB RAM). So for projects with normal complexity (cg=1.9) and even large projects with high complexity (cg=3) computation time of MCS increases linear by number of activities and complexity as long as activities are not correlated. So 2 seconds for a large simulation seems to be an acceptable time. If you like to study the effect of correlated activities the mean cpu - time is required for generation of a large amount of correlated random numbers (cholesky decomposition or something else), so that computation time exhibit quadratic growth. The analytical approximation is faster about a factor 10 and correlation is of no concern because the algorithm always processes a complete covariance matrix. The advantage of speed is diminished by higher complexity because the number of evaluations of moments in max – operations increase by a square law if you like to minimize skewness as proposed for the algorithm. Furthermore the algorithm exhibits quadratic growth due to the updates of covariance matrix.

Stefan Isermann / Critical Path Analysis

27.08.2015

41

5. Algorithms Conclusion We have seen that nowadays the response time of simple MCS may be sufficient provided it is carefully implemented so that it is in the range of a few seconds. The analytical approximation gives much more insight and is faster about a factor 10 depending on complexity of the project. Estimated precision of mean and uncertainty (standard deviation) of project duration is completely sufficient for project estimations and correlation of activities is incorporated. Although I have developed all necessary formulas to incorporate third order cumulants (s. appendix A.5), it is not advised to implement a suitable algorithm to a get a better precision for skewed distributions (1-2% percent for standard deviation!) at the expense of much more computation time which is required for the control of all third order cumulants! So you may choose between the normal approximation presented here or a simple MCS.

Stefan Isermann / Critical Path Analysis

27.08.2015

42

CONTENTS 1. Critical Path Analysis – Deterministic Methods 2. Uncertainty and probabilities and how they affect project duration 3. Main results from statistical analysis 4. Criticality and severity of activities 5. Algorithms A Mathematical Appendix References

Stefan Isermann/ Critical Path Analysis

27.08.2015

43

CONTENTS A A.1 A.2 A.3 A.4 A.5 A.6 A.7 A.8

Mathematical Appendix Distribution of the sum of random variables Covariance matrix of correlated paths Exact results for the maximum of two normal random variables Approximation error due to skewness Correction of moments of maximum due to skewness 2nd and 3rd order cumulants in critical paths analysis Number of possible paths in a project plan Simulation of dependencies with Copulas Stefan Isermann/ Critical Path Analysis

27.08.2015

44

Mathematical Appendix A.1 Distribution of the sum of random variables If you add two independent random variables with density f(x1) and g(x2) respectively the resulting density of the sum h(x=x1+x2) is given by the convolution of f and g (A1.1) h( x) = ( f ∗ g )( x) =

∞

∞

−∞

−∞

∫ f (t ) g ( x − t )dt =
∫ f ( x − t ) g (t )dt


In the case of dependent variables you have to replace f(x1) and g(x2) by a common density f(x1, x2) and you get h( x) =

∞

∫

−∞

∞

f (t , x − t )dt =
∫ f ( x − t , t )dt
−∞

(A1.2)

If the variables are normal distributed the density of the sum is again normal distributed with µ = µ 1 + µ 2 and σ 2 = σ12 + 2ρσ1σ2 + σ22. That the sum of gaussians is a gaussian once again is a feature of the gaussian distribution and does not hold generally for other distributions. But due to the central limit theorem the sum of any skewed random variables tend to a gaussian distribution, especially the skewness of the sum of n independent random variables shrinks by n1/2 (Berry –Esseen theorem). In A.5 we show that skewness may also be reduced by max – operations (excluding the pathological cases) , so that critical path distribution should always be nearly gaussian for n not too small. Stefan Isermann / Critical Path Analysis

27.08.2015

45

Mathematical Appendix A.2 Covariance of correlated paths Assume path x consists of n activities path y of m activities and they have r activities in common. The activities are independent and have standard deviation σi. The covariance of the path x and path y is   cov(x , y) = cov ∑ x , ∑ y  = (A2.1) ∑σ   For a project of n activities, which lead to k paths in the project, we define a path-activity n x k - matrix p with n

i =1

1 pi , j =  0

m

i

j=1

r

i

2 c c = common .path

if activity i is member of path j else

(A2.2) Each activity has a standard deviation σi which defines a diagonal n x n covariance matrix of activities σ i2 , i = j (A2.3) Σ activites =  0, i ≠ j

Then the k x k covariance matrix of project paths Σpaths is given by Σ paths = p T Σ activitiesp

(A2.4) where T means the transpose and Σjj,paths equals the variance of the jth path. Obviously the non-diagonal elements do not vanish if paths have activities in common.

Stefan Isermann / Critical Path Analysis

27.08.2015

46

Mathematical Appendix A.2 Covariance matrix of correlated paths If activities i and j are correlated among themselves with correlation ρij the n x n covariance matrix of activities is not diagonal Σ activites = ρijσi σ j (A2.5) and Σactivities in (A2.4) has to be replaced by (A2.5). ij

Important note: If you like to introduce correlations between activities (from where do you know??) via (A2.5) you have to be careful. You must not choose correlations as you like, because a symmetric matrix is not necessarily positive definite, but correlation matrix has to be. A Cholesky decomposition may work as a test. If it fails your correlation matrix is not positive definite and the resulting matrix in (A2.5) is not a covariance matrix.

Stefan Isermann / Critical Path Analysis

27.08.2015

47

Mathematical Appendix A.2 Covariance matrix of correlated paths A quick and easy way to create an appropriate correlation matrix is to choose a constant positive correlation between all activities. Another way is to choose different positive correlations but constant in a block. E. g. if ρ2,3 = r and ρ2,4 = r then ρ3,4 should also equal r (r ≥ 0). Figure A2.1 shows the restricted region of possible values of a 3x3 correlation matrix. The probability to chance on a valuable correlation matrix is about 62% (π2/16). In higher dimensions the region will further reduce, so it becomes less likely to get allowed values by random choice (the inner region is a projection in three dimensions). To generate random correlation matrices, you may take random numbers for the elements of a lower triangular matrix L and divide them by the Euclidean norm of their row. Then LLT is a correlation matrix. If you are interested in creating realistic correlation matrices you should read the paper of Hardin, Garcia [6].

Fig. A2.1

Stefan Isermann / Critical Path Analysis

27.08.2015

48

Mathematical Appendix A.2 Covariance matrix of correlated paths At least we like to give some general remarks about the impact of correlated activities on the critical path. If you introduce positive correlations between activities particularly the same correlation to all activities it will reduce the shift of the maximum and so the shift of the critical path! This is due to the fact that correlation between activities will broaden the path distribution but enlarge the normalized distance. Look at θ (eq. (3.4)) of two neighboring paths. In the uncorrelated case we get θ = ∑σ (A2.6) whereas with a correlation ρ for all activities you get θ = ∑σ − 2ρ ∑σ σ (A2.7) So θ is reduced by the covariance of the activities of the two paths which are not in common. Due to the positive correlation of the activities the critical path probability of the path with smaller mean decreases. Nevertheless you may construct relations where the shift increases by introducing some correlations. For example the exclusive introduction of correlations between activities which are not in common will broaden the path distribution but will not alter the covariance of the paths. This will lead to a smaller normalized distance but seems to be rather artificial. 2

2 i i∈{not common activities}

2

2 i i∈{not common activities}

i j i , j∈{ not common activities}

Stefan Isermann / Critical Path Analysis

27.08.2015

49

Mathematical Appendix A.3 Exact results for the maximum of two normal random variables You will find the exact results for the maximum of two gaussian random variables in several papers, for example [1], [7] and [8]. Whereas the paper of C.E. Clark [1] first presents these formulas to approximate the moments of the maximum of a finite set of random variables. With definitions θ = σ + σ − 2 ρσ σ (A3.1) µ −µ α= (A3.2) θ the probabilities P(x1> x2) and P(x2> x1) are (A3.3) P( x > x ) = Φ (α ) = 1 − P( x > x ) = 1 − Φ (−α ) and the moments of the maximum of two gaussians X1 and X2 with means µ1, µ2, standard deviations σ1, σ2 and linear correlation coefficient ρ are ν = µ Φ (α ) + µ Φ (−α ) + θϕ (α ) (A3.4) ν = (µ + σ )Φ (α ) + (µ + σ )Φ ( −α ) + (µ + µ )θϕ (α ) (A3.5) ν = (µ + 3µ σ )Φ (α ) + (µ + 3µ σ )Φ ( −α ) + (A3.6) 2

2

2

1

2

1

1

2

1

1

1

2

2

1

2

2

2

1

3

3 1

(

1

2

2

1

2

2

2

2

1

1

1

3

) (

2

2

2

2

2

(

2 4 2 2  2 2σ 1 + σ 1 σ 2 + 2σ 2 − 2 ρσ 1σ 2 σ 2 + σ 2  µ1 + µ1µ 2 + µ 2 2 +  θ2  4

2

))θϕ (α )  

Stefan Isermann / Critical Path Analysis

27.08.2015

50

Mathematical Appendix A.3 Exact results for the maximum of two normal random variables If another gaussian random variable Y is linear correlated to X1 and X2 with coefficient ρ1 and ρ2 respectively the covariance of Y and maximum(X1, X2) has to be updated according to Cov(Y , max( X 1 , X 2 ) = Cov(Y , X 1 )Φ (α ) + Cov(Y , X 2 )Φ (−α ) or

(A3.7)

ρσ yσ max = ρ1σ yσ 1Φ (α ) + ρ 2σ yσ 2Φ (−α )

(

σ max = ν 2 −ν 12

)

1/ 2

So the covariance of Y and the maximum X1 and X2 is given by the weighted sum of the individual covariances according to their probability to be the maximum. To be complete we give the probability density function of the maximum Z=max(X1, X2) which reduces to eq. 2.5 for ρ = 0 ϕ (z) =

1

σ1

 z − µ1   1 Φ 2   σ1   1 − ρ

ϕ 

(

)

1/ 2

 z − µ2 z − µ1   1  z − µ 2   1  Φ  + ϕ  −ρ   σ1  σ 2  σ 2   1− ρ 2  σ2

(

)

1/ 2

 z − µ1 z − µ 2     −ρ σ 2    σ1

Stefan Isermann / Critical Path Analysis

(A3.8)

27.08.2015

51

Mathematical Appendix A.4 Approximation error due to skewness Shina et. al. [7] introduce the error of the gaussian approximation of the pdf of the maximum of two normal distributed variables by (A4.1) Ξ = ϕ ( x ) − ϕ ( x) dx ∞

Z , ZG

∫

Z

ZG

−∞

Where ϕZG denotes the gaussian approximation and ϕZ the correct pdf of the maximum ( s. (3.4)). It turns out, that this error depends on three canonical parameters (normalized distance α, correlation ρ and ratio of standard deviations σ=σ1/σ2). They calculate this error numerically depending on these three parameters and stored the values up to any desired precision in a LookUp–table for minimizing the error by succeeding combination of suitable variables. This heuristic method works fine but is not necessary. Skewness is invariant under linear transformations and therefore the natural measure of the goodness of the gaussian approximation of the correct pdf of the maximum. Furthermore skewness is directly correlated to the error defined in eq. (4.1). To show this we do a cumulant expansion of the skewed normal distribution known as Edgeworth series taking into account the skewness γ. In one dimension we know that the first order of the Edgeworth expansion for pdf f (x) and Stefan Isermann / Critical Path Analysis

27.08.2015

52

Mathematical Appendix A.4 Approximation error due to skewness cdf F(x) is

1 1  1  f(x) = ϕ(x) − γ( 3 x − x 3 )ϕ(x),with:ϕ(x) = exp − x 2  6 2π  2 

(A4.2)

x

1 F(x) = Φ(x) − γ(x 2 − 1 )ϕ(x),with:Φ(x) = ∫ ϕ(t)dt 6 −∞

If we calculate the integral (4.1) with ϕZ (x) ≈ f(x) we find Ξ Z ,ZG ≈

4e − 3 2 + 1 γ = 0.252γ 3 2π

(A4.3) which means, that the approximation error defined in [7] is approximately γ/4 as long as the exact cdf is not heavily skewed. Furthermore skewness arises directly by calculating first and second moment of the Gaussian approximation (s. A.3), so you do not have to calculate any integral. The following figures show this relation directly. Fig. 4.1 shows a plot of skewness for ρ = 0.5 depending on α, σ. Fig. 4.2 shows contour lines of skewness for different values of ρ. Fig. 4.3 and 4.4 show correspondent plots for the numerically calculated error Ξ. Colors indicate the region where the approximation is fairly good (white to green) and where it fails (red). The structure of both is the same. For small γ you will find the same contour lines in the skewness plot as in the error plot deviding the contour level by approximately 4! If you calculate the error ΞZ,ZG numerically you have to take care of the singularity at σ = 0 where ρ is not defined and the pdf turns to a delta – distribution! Sinha et. al. do not seem to consider this fact as you can see from fig. 4 in [7] (you may compare this figure with fig. 4.3 here). Stefan Isermann / Critical Path Analysis

27.08.2015

53

Mathematical Appendix A.4 Approximation error due to skewness

γ1 Fig. A4.1: skewness (α,σ) for ρ=0.5

Fig. A4.3: error Ξ (α,σ) for ρ=0.5

Fig. A4.2: contour lines of skewness (α,σ) for different values of ρ

Fig. A4.4: contour lines of error Ξ (α,σ) for different values of ρ

54

Mathematical Appendix A.5 Correction of moments of maximum due to skewness As in the one dimensional case we consider the impact of skewness on the moments of the maximum with the help of the bivariate Edgeworth 1 ∂ (A5.1) f ( x , x ) = ϕ ( x , x ) + ∑ (−1) A ϕ(x , x ) ∞

1

2

1

2

n+ m

n+m

n + m ≥3

nm

n! m! ∂x1n ∂x2 m

1

2

where Anm are functions of the cumulants and ϕ(x1,x2) is the bivariate normal distribution. For simplicity we consider the impact of γ1, which means n=3, m=0 and A30=γ1σ13. We get f ( x1 , x2 ) ≅

1 2πσ 1σ 2

(

x − µ1 ~ x2 − µ 2 , x2 = with : ~ x1 = 1

σ1

)(

 1 ~ exp − x12 − 2 ρ~ x1 ~ x2 + ~ x22 2 2 2 1 − ρ 1− ρ 

  1 1 1 − γ 1  6 1 − ρ 2 

)

(

2 ~   ~ ~   ~  x1 − ρx2  3 −  x1 − ρx2    3/ 2  1 − ρ 2   1 − ρ 2        

)

(A5.2)

σ2

For the uncorrelated case you easily identify the Edgeworth expansion (4.2) in x1 multiplied by the normal distribution of x2 which indicates the desired marginal distributions of f(x1,x2).

Stefan Isermann / Critical Path Analysis

27.08.2015

55

Mathematical Appendix A.5 Correction of moments of maximum due to skewness The correction of the cdf of the maximum of x1, x2 due to skewness is z z

Fγ 1 ( z = max( x1 , x2 )) =

∫ ∫ fγ (x , x )dx dx 1

1

2

1

2

− ∞− ∞

 z − ρz   z12 − 2 ρz1 z 2 + z 22    1  z − ρz 2  1 1   −   = − γ 1  z12 − 1 ϕ ( z1 )Φ 2 + ρ  z1 + 1 exp  1− ρ 2  6  1 − ρ 2  2π 1 − ρ 2 2 1− ρ 2        z − µ1 z − µ2 with : z1 = , z2 =

(

)

σ1

(

)

(A5.3)

σ2

Once more for ρ = 0 you may identify the cdf of the correspondent Edgeworth Expansion of z1 (s. (4.2)) multiplied by the cdf of z2. Furthermore you may easiliy verify that Fγ ( z = ∞) = 0 (the total probability has to equal 1) so that P (x > x ) = −P (x > x ) (A5.4) and 1

γ1

2

γ1

1

Pγ 1 ( x2 > x1 ) =

1

2

∞ x2

∫ ∫ fγ (x , x )dx dx 1

1

2

1

2

− ∞− ∞

1 σ  = − γ 1  1  α 2 − 1 ϕ (α ) 6 θ  3

(

)

with : θ = σ 12 − 2 ρσ 1σ 2 + σ 22 , α =

(A5.5) µ1 − µ 2 θ Stefan Isermann / Critical Path Analysis

27.08.2015

56

Mathematical Appendix A.5 Correction of moments of maximum due to skewness So if α < 1, that is when paths are neighbored and skewness becomes larger by combining their normal approximations, the skewness of the probability distribution of a path will reduce its probability to be the maximum compared to the normal approximation and enhance the probability of the adjacent path! By interchanging the indices 1 and 2 we get 1  σ  σ  Pγ ( x2 > x1 ) =  γ 2  2  − γ 1  1   6  θ  θ  3

3

 2  α − 1 ϕ (α )  

(

)

(A5.6)

Pγ ( x1 > x2 ) = − Pγ ( x2 > x1 )

It is bit tedious to calculate the correction to the momentum generating function Mγ1(t) ∞

M γ1 ( t ) = ∫ e tz −∞

∂Fγ1 ( z ) ∂z

(A5.7)

∞

dz = − t ∫ e tz Fγ1 ( z )dz −∞

From where you may get the correction of the moments νn induced by skewness γ1 via ν n γ1 =

∂ n M γ1 (t = 0) ∂t n

=n

∂ n −1m γ1 (t = 0) ∂t n −1

∞

, with : m γ1 ( t ) = − ∫ e tz Fγ1 ( z )dz −∞

(A5.8) Stefan Isermann / Critical Path Analysis

27.08.2015

57

Mathematical Appendix A.5 Correction of moments of maximum due to skewness So we have to calculate the integral using eq. (A5.3) ∞

mγ 1 (t ) = − ∫ etz Fγ 1 ( z )dz −∞

3 3 γ1   σ1   σ1  2 =  −   ( µ1 − µ2 ) +   σ1 − 3 ρσ1σ 2 + 2σ 22 6   θ  θ

(

+

)t ϕ(α ) exp µ σ (σ  

1 2



2

(

)

− ρσ1 ) + µ2 σ1 (σ1 − ρσ 2 ) 1 − ρ 2 σ12 σ 22 2  t+ t  2 θ 2θ 2 

γ1  3 2  σ (σ − ρσ 2 )  1   σ1 t Φ  α + 1 1 t  exp µ1t + σ12t 2   6 θ 2    

(A5.9)

Introducing the following abbreviations a1 = µ1 − µ2 , b1 = σ12 − 3 ρσ1σ 2 + 2σ 22 ,

µ1σ 2 (σ 2 − ρσ1 ) + µ2 σ1 (σ1 − ρσ 2 ) , θ2 1 − ρ 2 σ12 σ 22 d= 2θ 2

c=

(

(A5.10)

)

Stefan Isermann / Critical Path Analysis

27.08.2015

58

Mathematical Appendix A.5 Correction of moments of maximum due to skewness we get with (A5.8) the correction to the expected mean of the maximum 3

ν1γ1 = −

γ1  σ1    ( µ1 − µ2 )ϕ (α ) 6θ 

(A5.11)

The correction for the second moment γ 1  σ1    (b1 − ca1 )ϕ (α ) 3θ  3

ν 2γ = 1

(A5.12)

And the correction for the third moment 3

ν3 γ1 =

(

(

) )

γ1  σ1  2 3   2b1c − c + d a1 ϕ (α ) + γ1σ1 Φ (α ) 2θ 

(A5.13) The corrections due to γ2 are given by interchanging the indices 1 and 2 and replacing α by –α.

Stefan Isermann / Critical Path Analysis

27.08.2015

59

Mathematical Appendix A.5 Correction of moments of maximum due to skewness Fig. A5.1 shows the effect of a skewed compared to a normal distribution. The upper left shows the bivariate 1. order Edgeworth Expansion (A5.2) with µ1=0, σ1=0.8, γ1=0.6, µ2=0.2, σ2=0.9, γ2=0 and ρ=0.5. The upper right is the bivariate normal distribution with same parameters but γ1=0. The lower left and right shows the corresponding maximum probability distributions (green) and cumulative distributions (red) and in dashed lines the conditional probabilities of x1 and x2 respectively. The moments µ, σ, γ of the maximum and the probability p1for x1 > x2 (filled area) are obtained via numerical integration and exactly coincide with the values given by eqs. A3.1 – 3.6 and A5.6, A5.11 – 5.13. It is obvious that correction due to skewness becomes major for higher moments.

Fig. A5.1: Edgeworth Expansion

Stefan Isermann / Critical Path Analysis

27.08.2015

60

Mathematical Appendix A.5 Correction of moments of maximum due to skewness To see the effect of the Edgeworth expansion, we compare the different approximations to the maximum distribution of two independent and skewed beta distributed random variables (red and blue). The exact cumulants µ, σ and γ of the maximum distribution (green) and probability p1 (red > blue) in the upper figure are obtained via numerical integration of the known pdf. The normal approximation neglecting the skewness of the beta distributions gives reasonable approximations of mean and variance, especially the mean is nearly exact. Values are obtained via formulas A3.3 - 3.5. Taking skewness of the betas into account gives nearly an exact value for variance too whereas skewness slightly differs from the exact value. Values are obtained via A3.3 - A3.6 combined with A5.6, A5.11 A5.13 (beta dist.: a=2, b=9, α=1.5, β=4.5, a=3, b=7, α=2, β=4). For most practical reasons where only an estimation of variance of activities is done, the normal approximation should be sufficient.

Fig. A5.2: comparison of approximations

Stefan Isermann / Critical Path Analysis

27.08.2015

61

Mathematical Appendix A.5 Correction of moments of maximum due to skewness Up to here we have neglected third order cumulants κ2,1=γ21σ12σ2 and κ1,2=γ12σ1σ22. But as activities become member of different paths, the paths will interact and generate third order cumulants in the same way as they become correlated (s. appendix A.2 and A.6). Taking into account third order cumulants you have to do a lot of integrals or derivatives, because f ( x1 , x2 ) = ϕ ( x1 , x2 ) +

∞

∑ (−1) n+ m Anm

n + m ≥3

∞  1 ∂ n+m 1 ∂ n+m   ϕ ( x1 , x2 ) ϕ ( x , x ) = 1 + A nm 1 2  n∑ n!m! ∂x1n ∂x2 m n!m! ∂µ1n ∂µ 2 m  + m ≥3 

(A5.14)

Here are the results ν1γ21 =

1 σ 12σ 2 γ21 3 ( µ1 − µ2 )ϕ (α ) 2 θ

ν2 γ21 = γ21 ν3 γ21 =

with

σ 12σ 2 θ3

(ca1 − b21 )ϕ (α )

(A5.15)

3 σ 12σ 2 2 γ21 3 c + d a1 − 2b21c ϕ (α ) 2 θ

((

)

)

b21 = σ 22 − ρσ 1σ 2

and a1, c and d as in eqs. (A5.10) The corrections due to γ12 are given by interchanging the indices 1 and 2, whereas γijσ1iσ2j are the values of the third order cumulants of X1 and X2. Stefan Isermann / Critical Path Analysis

27.08.2015

62

Mathematical Appendix A.5 Correction of moments of maximum due to skewness It remains to compute the effect of κ2,1 and κ1,2=γ12σ1σ22 on the maximum probabilities Pγ 21 ( x2 > x1 ) = Pγ 21 ( x1 > x2 ) =

∞ x2

1 σ 12σ 2 2 f ( x , x ) dx dx = − γ ∫−∞−∫∞ γ 21 1 2 1 2 2 21 θ 3 1 − α ϕ (α )

(

∞ x1

)

(A5.16)

1 σ 12σ 2 2 f ( x , x ) dx dx = γ ∫−∞−∫∞ γ 21 1 2 2 1 2 21 θ 3 1 − α ϕ (α )

(

)

So we get for the probabilities 1 σ3 σ3 σ 2σ σ σ2  P( x1 > x2 ) = φ (α ) +  γ 2 32 − γ 1 13 + 3γ 21 1 3 2 − 3γ 12 1 2  1 − α 2 ϕ (α ) θ θ θ  6 θ P( x2 > x1 ) = 1 − P( x1 > x2 )

(

)

(A5.17)

And at least the third order cumulants involving the maximum of two gaussian variables κ (t , max( x1 , x2 ), max( x1 , x2 )) =< t max( x1 , x2 ) 2 > −2(cov(t , x1 )φ (α ) + cov(t , x2 )φ (− α ))(µ1φ (α ) + µ 2φ (− α ) + θϕ (α ))

(

)ϕθ(α )

< t max( x1 , x2 ) 2 >= 2µ1 cov(t , x1 )φ (α ) + 2µ 2 cov(t , x2 )φ (− α ) + 2 cov(t , x1 )σ 12 + cov(t , x2 )σ 22 − cov( x1 , x2 )(cov(t , x1 ) + cov(t , x2 ) )

and κ (t , t , max( x1 , x2 )) =

(cov(t , x1 ) − cov(t , x2 ))2 ϕ (α )

(A5.18) (A5.19)

θ

Stefan Isermann / Critical Path Analysis

27.08.2015

63

Mathematical Appendix A.6 2nd and 3rd order cumulants in critical paths analysis Let us assume we add duration of successor C to duration of activity A and activity B. Then we get from multilinearity of cumulants cov( A + C , B + C ) = cov( A, B + C ) + cov(C , B + C ) = cov( A, B) + cov( A, C ) + cov( B, C ) + cov(C , C ) (A6.1) Even if all activities are uncorrelated new random variables A’ = A + C and B’ = B + C becomes correlated with cov(A’,B’) = var(C). So we have to consider the general case of correlated random variables, where in each step of adding durations the covariance matrix has to be updated according to cov( X i + X k , X j ) = cov( X i , X j ) + cov( X k , X j ) cov'ij = cov ij + cov kj , i ≠ j , k ≠ i

(A6.2) And in the case where Xk is the maximum of n different normal random variables Xr which occur with probability pr n

cov( X i + X max , X j ) = cov( X i , X j ) + ∑ pr cov( X r , X j ) r n

cov'ij = cov ij + ∑ pr cov rj , r ≠ i, i ≠ j

(A6.3)

r

Stefan Isermann / Critical Path Analysis

27.08.2015

64

Mathematical Appendix A.6 2nd and 3rd order cumulants in critical paths analysis The diagonal elements of the updated covariance matrix are n

cov'ii = var( X i + X max ) = var( X i ) + var( X max ) + 2∑ pr cov( X r , X i )

(A6.4) If you take skewness into account you have to look at all third order cumulants additionally r

κ ( X i + X k , X i + X k , X i + X k ) = κ ( X i , X i , X i ) + 3κ ( X i , X k , X k ) + 3κ ( X i , X i , X k ) + κ ( X k , X k , X k ) κ (Xi + Xk , X j , X j ) = κ(Xi, X j , X j ) +κ (X k , X j , X j )

(A6.5)

κ ( X i + X k , X i + X k , X j ) = κ ( X i , X i , X j ) + κ ( X k , X k , X j ) + 2κ ( X i , X k , X j )

where you may neglect 3-activity interaction (κ (Xi, Xj, Xk) = 0) in the last line. Nevertheless if you do an Egdeworth expansion you have to compute third order cumulants containing the maximum of two random variables κ (Xi, Xi, max( Xj, Xk)) and κ (Xi, max( Xj, Xk) , max( Xj, Xk)) , which is a bit tedious (see results in previous chapter).

Stefan Isermann / Critical Path Analysis

27.08.2015

65

Mathematical Appendix A.7 Number of possible paths in a project plan As we have seen in chapter 1 we get the number of paths at a node by just adding all paths of the predecessors. This will directly lead to the generalized theory of Fibonacci numbers. First let us define the complexity cg of the graph as number of all edges c = (A7.1) number of nodes - 1 where the number of nodes include the starting and the end point. So cg is just the average number of edges leaving a node. For the project example in chap. 1 we have cg = 16/(11-1) = 1.6. Obviously the number of paths should be influenced by the complexity in the same way as the number of paths at each node increases due to the number of predecessors. That’s the reason why we take a closer look at Fibonacci numbers defined by Fn+2= Fn+1+ Fn. We can write this recursive formula as a matrix equation  F  1 1  F  1 1   F  (A7.2)   =   =      g

n +1

n+2

n +1

 Fn +1  1 0  Fn  1 0 

1

 F0 

or F  1 1   , Fk =  k +1  Fn = A n F0 with A =  1 0   Fk 

(A7.3) Stefan Isermann / Critical Path Analysis

27.08.2015

66

Mathematical Appendix A.7 Number of possible paths in a project plan So Fn will increase for large n due to the largest eigenvalue of A, which we get from the characteristic equation of A (λ − 1)λ − 1 = 0 λ1 =

1+ 5 1− 5 , λ2 = 2 2

(A7.4) which means that in a simple project structure, where each node has two direct predecessors (in a topological sorted DAG) the number of paths increases proportional to 1.618 n. And this is the maximum rate at which the number of path in a project of complexity cg= 2 may increase. In a more general formulation , the recursion F = ∑ c F , 0 ≤ c ≤ 1, ∑ c = c (A7.5) leads to a characteristic equation λ − ∑c λ = 0 (A7.6) which has always a positive real root between 1 and 2 (follows from Descartes rule of sign and Rouché theorem). So 2n is obviously the maximum rate the number of paths may increase, which is the case when every node in a topological sorted DAG has the maximum number of predecessors. Let us take a look at a more realistic graph which represent typical projects and which we use in our simulations . m

m +1

m

i =1

m

i

i

i =1

i

graph

i −1

m

i =1

i

i

Stefan Isermann / Critical Path Analysis

27.08.2015

67

Mathematical Appendix A.7 Number of possible paths in a project plan We generate a random topological sorted DAG by the way that each node j has i successors with probability pi where the i successors are taken from a random sample of j+1, …, j+n. If it happens that a node j has no predecessor we will add a random predecessor in the range j-1, …, j-n So each node j will have n predecessors with c edges at average

λ5 − 0.38 ⋅ (λ4 + λ3 + λ2 + λ1 + 1) = 0

we get an upper bound for the increasing rate with λmax = 1.26. The corresponding simulation shows that for projects with a large number N of activities, the number of paths may exceed the number oft atoms in our observable universe ( O(1080))!

1e+39 1e+18

For example p1=0.5, p2=0.3, p3=0.1, p4=0.1 und n=5 lead to c ≈ 0.38 with complexity of the graph cg= n⋅ c ≈ 1.9. From

1e-03

k ≤n 1 (1 − P )n , P = 1 ∑ i ⋅ pi n n i =1

1e+60

1e+81

ratereg = 1.25

number of paths

c= P+

project with complexity c g = 1.9 , range n = 5

0

200

400

600

800

1000

N

Fig. A7.1: number of paths in a project

Stefan Isermann / Critical Path Analysis

27.08.2015

68

Mathematical Appendix A.8 Simulation of dependence with Copulas If you like to do MCS with correlated and not normal distributed activities than you have to consider copulas. From Sklar’s theorem it follows that for contious distributions, the univariate margins and the multivariate dependence structure can be totally separated, with the dependence structure fully explained by a copula (couples the dependence structure with marginal distributions). Sklar’s Theorem: Let H be an n-dimensional distribution function with margins F1,…,Fn. Then there exists an n-copula C such that for all x in Rn H ( x ,...., x ) = C (F ( x ),..., F ( x ) ) (A8.1) The copula is unique if margins F1,…,Fn are all continuous. So for any u in [0,1]n and continuous margins C (u ,...., u ) = H (F (u ),..., F (u ) ) (A8.2) As an example we give the Gauss – Copula with correlation ρ and standard bivariate normal distribution function Φ2 C (u , u ) = Φ (Φ (u ),Φ (u ), ρ ) (A8.3) and the bivariate copula density ϕ (Φ (u ),Φ (u ), ρ ) (A8.4) c (u , u ) = ϕ (Φ (u ) )ϕ (Φ (u ) ) 1

2

1

1

−1 1

2

1

2

2

ρ

2

1

2

n

−1 n

n

−1

1

−1

Ga

n

1

−1

Ga

ρ

1

2

−1

1

−1

2

−1

1

2

Stefan Isermann / Critical Path Analysis

27.08.2015

69

Mathematical Appendix A.8 Simulation of dependencies with Copulas Another example of an elliptical copula is the tν - copula with ν degrees of freedom and correlation parameter ρ C (u , u ) = t (t (u ), t (u ),ν , ρ ) (A8.5) -1 where t2 is the bivariate t-distribution and t the inverse of the univariate t–distribution withν degrees of freedom. The figure shows the densities of a Gauss - and a t4 - copula with parameter ρ = 0.5. The Gauss – copula exhibits no tail dependency, that is no correlation of extreme events. t

ρ ,ν

1

1

−1 2 ν

−1

1

ν

2

It is known that the Pearson’s (linear) correlation coefficient ρP is only invariant under linear transformation whereas the rank correlations Spearman’s ρS and Kendall’s τK are invariant under strictly monotonic transformations and are an intrinsic feature of the copula. Fig. A8.1: elliptical copulas Stefan Isermann / Critical Path Analysis

27.08.2015

70

Mathematical Appendix A.7 Simulation of dependencies with Copulas For Gaussian copula it is known that Spearman’s ρS is ρs =

ρ  arcsin P  π  2  6

where ρP is the parameter of the Gaussian copula in eq. A8.3. And generally it holds that 1 1

ρ s = 12∫ ∫ u1u2 dC (u1 , u2 ) − 3 = 0 0

< u1u2 > −1 / 4 Cov (u1u2 ) = = ρ P (u1 , u2 ) 1 / 12 Var (u1 ) Var (u1 )

(A8.6) (A8.7)

which means that linear correlation of copula variables u1 and u2 equals Spearman’s correlation ρS. Furthermore for copulas from elliptical distributions like Gauss - and tν – distribution, Kendall’s τK is also fully determined by the copula parameter ρP 2 τ = arcsin(ρ ) (A8.9) π For ρP = 0.5 we get for the Gauss copula ρS = 0.483 and Gauss and tν - copula τK = 1/3. With the density of the dependence structure generating copula the density of a multivariate distribution is given by the copula density multiplied with the marginal densities. For example a density from a Gauss copula: ϕ (Φ (F ( x ) ),Φ (F ( x ) ), ρ ) h( x , x ) = f (x ) f (x ) ϕ (Φ (F ( x ) ))ϕ (Φ (F ( x ) )) (A8.10) K

P

−1

2

1

2

−1

1

1

−1

2

2

−1

1

1

1

2

1

2

2

2

Stefan Isermann / Critical Path Analysis

27.08.2015

71

Mathematical Appendix A.8 Simulation of dependencies with Copulas Fig. A8.2 shows two bivariate beta distributions with marginal distributions β1(3,3) and β2(2,4). The structure is different depending on their copula and has a linear correlation ρP different from parameter ρ = 0.5 of the generating copula. Simulating from a Gauss (or tν - copula) is simple but cpu time consuming (see fig. A8.3). 1. Do a Cholesky decomposition C of correlation matrix R 2. Generate n random independent zi from N(0,1)) 3. Set y=Cz 4. Set ui=Φ(yi) 5. Set xi=Fi-1(ui) For more details see reference [3]. Fig. A8.2: Bivariate β distributions from different copulas with identical margins Stefan Isermann / Critical Path Analysis

27.08.2015

72

Mathematical Appendix A.8 Simulation of dependencies with Copulas Fig. A8.2 shows different amount of computation time necessary for simulation of activities and dependent of marginal distributions and correlations. You see that simulating with Copulas is rather time consuming.

Fig. A8.3 Stefan Isermann / Critical Path Analysis

27.08.2015

73

References [1] C. E. Clark, The greatest of a finite set of random variables, in Operations Research, Vol. 9, No. 2 (Mar – Apr 1961), pp. 145-162 [2] R. B. Arellano-Valle, M.G. Genton, On the exact distribution of the maximum of absolutely continuous dependent random variables, in Statistics & Probability Letters, Vol. 78, No. 1 (2008), pp. 27-35 [3] P. Embrechts, F. Lindskog, A. McNeil, Modelling Dependence with Copulas and Applications to Risk Management, in Handbook of Heavy Tailed Distributions in Finance, ed. S. Rachev, Elsevier, Chapter 8, pp. 329-384 [4] A. Genz, Comparison of Methods for the Computation of Multivariate Normal Probabilities , in Computing Science and Statistics Vol. 25 (1993), pp. 400-405 [5] A. Genz, Numerical Computation of Multivariate Normal Probabilities, in J. Comp. Graph Stat. 1 (1992), pp. 141149 [6] J. Hardin, S. R. Garcia, D. Golan, A method for generating realistic correlation matrices, in The Annals of Applied Statistics, Vol. 7, No. 3 (2013), pp. 1733 -1763 [7] S. Nadarajah, S. Kotz, Exact distribution of the max/min of two Gaussian random variables, in IEEE Transactions on Very Large Scale Integration (VLSI) Systems, Vol. 16, No. 2 (Feb 2008), pp. 210-212

Stefan Isermann / Critical Path Analysis

27.08.2015

74

References [8] D. Sinha, H. Zhou, and N. Shenoy, Advances in Computation of the Maximum of a Set of Gaussian Random Variables, IEEE Transactions on Computer-Aided Design, 26(8) (Aug 2007), pp. 1522-1533 [9] M.-J. Yao, W.-M. Chu, A new approximation algorithm for obtaining the probability distribution function for project completion time, Computers and mathematics with Applications, 54 (2007), pp. 282-295 [10] S. Tattoni, M. Schiraldi, L. Laura, Estimating project duration in uncertain environments: Monte Carlo Simulation strikes back, 22nd IPMA World Congress “Project Management to Run”, 9 - 11 November 2008, Roma, Italy [11] R Core Team (2015), R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. http://www.R-project.org/

Stefan Isermann / Critical Path Analysis

27.08.2015

75