print(`PROCEDURES:`): print(`Queue(X, Y, n, k, s), Cov(X, Y, a, b), kCov(X, Y, a, b, n, k)`): print(`X is the distribution of time between arrivals.`): print(`Y is the service time distribution.`): print(`n is the total number of customers in the system.`): print(`k is the number of customers in the system at time 0.`): print(`s is the number of identical parallel servers.`): print(`a is the first customer of interest.`): print(`b is the second customer of interest (a < b).`): with(LinearAlgebra): # Queue(X, Y, n, k, s) # ----------------------------------------------------------------------------- # Computes the sojourn time distribution for the nth customer in an # M/M/s queue, given k customers are in the system at time 0. Queue calls the # subprocedure MMsQprob(n, k, s) (and subsequently calls Q(n, i, k, s)) # which recursively calculates the required probabilities of the nth # customer seeing exactly i customers, including himself, in the system upon # arrival for i = 1 to n + k customers. Calculations are based on # algorithms provided in Kelton and Law (1985). Calls the subprocedure # BuildDist(X, Y, n, k, s) which builds the conditional sojourn time # distribution for i = 1 customer to i = n + k customers in the system. # Queue mixes the probabilities with the conditional sojourn time # distributions to return the exact PDF of the sojourn time in the # APPL "list of lists format." Requires the subprocedures mentioned # above along with the APPL software. The exponential arrival and # service random variables must be defined in the APPL format. # The procedure call is Queue(X, Y, n, k, s), where X is the arrival # time distribution, Y is the service time distribution, n is the # customer of interest, k is the number of customers in the system at # time 0, and s is the number of identical parallel servers. Both X # and Y must be exponential random variables in the APPL list-of-lists # format. # Name : Queue.mw # Author : Billy Kaczynski # Language : MAPLE 9 # Latest Revision : 09/11/08 # ---------------------------------------------------------------------------- Queue := proc(X, Y, n, k, s) global rho; local i :: integer, lst :: list, TIS :: list; rho := 1 / Mean(X) / (s * 1 / Mean(Y)); lst := BuildDist(X, Y, n, k, s); TIS := Mixture(MMsQprob(n, k, s), lst); return TIS; end: # MMSQprob(n, k, s) # ------------------------------------------------------------------------------- # Computes Pk(n, i)'s for an M/M/s queue, which is the probability that # customer n will see i customers in the system counting himself at time # Tn with k customers initially in the system at time 0. Calls the subprocedure # Q(n, i, k, s) which recursively calculates the required probabilities using # the algorithms provided in Kelton and Law (1985). The procedure returns the # ordered probabilities for i = 1 customer (an empty queue) to i = n + k # customers in a list. Note that the parameter rho for an M/M/s queue is # rho = lambda / (s * mu). # Name : MMsQprob.mw # Author : Billy Kaczynski # Language : MAPLE 9 # Latest Revision : 09/03/08 # ------------------------------------------------------------------------------- MMsQprob := proc(n, k, s) local i :: integer, lst :: list; lst := []; for i from 1 to n + k do lst := [op(lst), Q(n, i, k, s)]; od; return lst; end: # Q(n, i, k, s) # ------------------------------------------------------------------------------- # Computes the single probability Pk(n, i) for an M/M/s queue recursively # according to the algorithms provided in Kelton and Law (1985). # Name : Q.mw # Author : Billy Kaczynski # Language : MAPLE 9 # Latest Revision : 09/03/08 # ------------------------------------------------------------------------------- Q := proc(n, i, k, s) option remember; global rho; local p, j :: integer, h :: integer; if (k >= 1) and (i = k + n) then if k >= s then p := (rho / (rho + 1)) ^ n elif k + n <= s then p := rho ^ n / (mul(rho + (k + j - 1) / s, j = 1 .. n)) elif (k < s) and (s < k + n) then p := rho ^ n / ((rho + 1) ^ (n - s + k) * (mul(rho + (k + j - 1) / s, j = 1 .. s - k))) fi; fi; if (k = 0) and (i = n) then if n <= s then p := rho ^ n / mul(rho + (j - 1) / s, j = 1 .. n) elif n > s then p := rho ^ n / ((rho + 1) ^ (n - s) * mul(rho + (j - 1) / s, j = 1 .. s)) fi; fi; if i = 1 then p := 1 - add(Q(n, j, k, s), j = 2 .. n + k) fi; if (k >= 1) and (i >= 2) and (i <= k) and (n = 1) then if k <= s then p := rho / (rho + (i - 1) / s) * mul(1 - rho / (rho + (k - j + 1) / s), j = 1 .. k - i + 1) elif (k > s) and (i > s) then p := rho / (rho + 1) ^ (k - i + 2) elif (i <= s) and (s < k) then p := rho / ((rho + 1) ^ (k - s + 1) * (rho + (i - 1) / s)) * mul(1 - rho / (rho + (s - j) / s), j = 1 .. s - i) fi; fi; if (n >= 2) and (i >= 2) and (i <= k + n - 1) then if i > s then p := rho / (rho + 1) * add((1 / (rho + 1) ^ (j - i + 1) * Q(n - 1, j, k, s)), j = i - 1..k + n - 1) elif i <= s then p := rho / (rho + (i - 1) / s) * (add((mul(1 - rho / (rho + (j - h + 1) / s), h = 1..j - i + 1)) * Q(n - 1, j, k, s), j = i - 1..s - 1) + product(1 - rho / (rho + (s - h) / s), h = 1..s - i) * add((1 / (rho + 1)) ^ (j - s + 1) * Q(n - 1, j, k, s), j = s .. k + n - 1)); fi; fi; return p; end: # BuildDist(X, Y, n, k, s) # ------------------------------------------------------------------------------- # Creates the appropriate conditional sojourn time distribution for each case # where a customer arrives to find i = 1 to i = n + k customers present in an # M/M/s queue with k customers intially present. The procedure call is # BuildDist(X, Y, n, k, s), where X is the arrival time distribution, Y is # the service time distribution, n is the customer of interest, k is the # number of customers in the system at time 0, and s is the number of # identical parallel servers. Both X and Y must be exponential random # variables in the APPL list-of-lists format. # Name : BuildDist.mw # Author : Billy Kaczynski # Language : MAPLE 9 # Latest Revision : 09/03/08 # ------------------------------------------------------------------------------- BuildDist := proc(X, Y, n, k, s) local i :: integer, lst :: list; lst := []; for i from 1 to n + k do if s = 1 then lst := [op(lst), ErlangRV(1 / Mean(Y), i)] else if i <= s or s > n + k then lst := [op(lst), Y]; else lst := [op(lst), Convolution(ErlangRV(s * 1 / Mean(Y), i - s), Y)]; fi; fi; od; return lst; end: #cases(n) #-------------------------------------------------------------------------------- # Generates all possible arrival/departure sequences for n customers in an M/M/1 # queue initially empty and idle. Resulting list of sequences consists # of 1's and -1's, where a 1 is an arrival and a -1 is a departure. # The sequences are returned in the matrix C, of dimension c by 2n, where # c in the nth Catalan number calculated as (2n)! / n! / (n + 1)!. The # procedure calls ini(n) to initialize the first sequence in the matrix, then # uses the procedures swapa(n, A), swapb(n, A), and okay(n, A) to # create the remaining sequences according to a prefix shift algorithm. # For each row in the resulting matrix, an associated path matrix can # be generated via the procedure path(n, A). # # Name : cases.mw # Author : Billy Kaczynski # Language : MAPLE 9 # Latest Revision : 01/12/09 #---------------------------------------------------------------------------------- cases := proc(n) options remember; local c, C, i; c := (2 * n)! / n! / (n + 1)!; C := Matrix(c, 2 * n); for i from 1 to c do if (i = 1) then C[[i], 1 .. -1] := ini(n); else C[[i], 1 .. -1] := swapa(n, C[[i - 1], 1 .. -1]); fi: if (okay(n, C[[i], 1 .. -1]) = 0) then C[[i], 1 .. -1] := swapb(n, C[[i - 1], 1 .. -1]); fi: od: return C; end: # kcases(n, k) # -------------------------------------------------------------------------------- # Generates all possible arrival/departure sequences for n customers in an M/M/1 # queue with k customers initially present. Resulting list of sequences consists # of 1's and -1's, where a 1 is an arrival and a -1 is a departure. # The sequences are returned in the matrix C. The procedure calls cases(n + k) # which susequently calls ini(n + k) to initialize the first sequence in the # matrix, then uses the procedures swapa(n, A), swapb(n, A), and okay(n, A) to # create the remaining sequences according to a prefix shift algorithm. C is # then simplified by deleting the rows where the first k entries are not 1's. # Furthermore, the first k columns are also deleted since they must all contain # 1's to represent the arrivals of the k customers present at time zero. # For each row in the resulting matrix, an associated path matrix can # be generated via the procedure kpath(n, k, A). # # Name : kcases.mw # Author : Billy Kaczynski # Language : MAPLE 9 # Latest Revision : 03/07/09 #---------------------------------------------------------------------------------- kcases := proc(n, k) local i, j, C; C := cases(n + k); j := 1; while j < RowDimension(C) + 1 do if (add(C[j, i], i = 1..k) <> k) then C := DeleteRow(C, j); else j := j + 1; fi: od: C := DeleteColumn(C, 1..k); return C; end: #ini(n) #-------------------------------------------------------------------------------- # Initializes the matrix C according to Ruskey and Williams (2008). Returns # the first row of C to enable use of their prefix shift algorithm. Requires # the parameter n, the number of customers. # # Name : ini.mw # Author : Billy Kaczynski # Language : MAPLE 9 # Latest Revision : 01/11/09 # ---------------------------------------------------------------------------------- ini := proc(n) options remember; local L, i; L := Matrix(1, 2 * n, -1): L[1, 1] := 1: for i from 3 to (n + 1) do L[1, i] := 1 od: for i from (n + 2) to (2 * n) do L[1, i] := -1 od: return L: end: # swapa(n, A) #-------------------------------------------------------------------------------- # Conducts the (k + 1)st prefix shift in creating all instances of the case # matrix, C, according to Ruskey and Williams (2008). Requires the # arguments n, number of customers, and A, row i of C. Returns the successor # of C[i, ] to be checked by the procedure okay(n, A). # # Name : swapa.mw # Author : Billy Kaczynski # Language : MAPLE 9 # Latest Revision : 01/11/09 # ---------------------------------------------------------------------------------- swapa := proc(n, A) local check, i, temp1, j, R; R := A; check := 1; for i from 2 to (2 * n - 1) do if ((R[1, i] = -1) and (R[1, (i + 1)] = 1)) then temp1 := R[1, i + 2]; R[1 .. 1, 3 .. (i+2)] := R[1 .. 1, 2 .. (i + 1)]; check := 0; R[1, 2] := temp1; fi: if (check = 0) then break fi; od: return R: end: # swapb(n, A) #-------------------------------------------------------------------------------- # Conducts the (k)th prefix shift in creating all instances of the case # matrix, C, according to Ruskey and Williams (2008). Requires the # arguments n, number of customers, and A, row i of C. Returns the successor # of C[i, ] in the event that okay(n, a) identifies the output of swapa(n, a) # invalid. # # Name : swapb.mw # Author : Billy Kaczynski # Language : MAPLE 9 # Latest Revision : 01/11/09 # ---------------------------------------------------------------------------------- swapb := proc(n, B) local check, i, temp, j; check := 1; for i from 2 to (2 * n - 2) do if ((B[1, i] = -1) and (B[1, (i + 1)] = 1)) then temp := B[1, i + 1]; B[[1 .. 1], [3 .. (i+1)]] := B[1 .. 1, 2 .. (i)]; check := 0; B[1, 2] := temp; fi: if (check = 0) then break fi; od: return B: end: # okay(n, E) #-------------------------------------------------------------------------------- # Checks the output of swapa(n, A) for an illegal prefix shift, meaning the # result contains an impossible arrival/service sequence. Requires arguments n, # number of customers, and E, resulting vector from swapa(n, A). If the # (k + 1)st shift is legal, okay(n, E) returns 1, signifying the correct # successor in C. If the (k + 1)st shift is illegal, okay(n, E) returns 0 # which in turn calls swapb(n, A) to produce the correct successor row in C. # # Name : okay.mw # Author : Billy Kaczynski # Language : MAPLE 9 # Latest Revision : 01/11/09 # ---------------------------------------------------------------------------------- okay := proc(n, E) local s, i, test; test := 1: s := 0; for i from 1 to (2 * n - 1) do s := s + E[1, i]; if (s < 0) then test := 0; break; fi: od: return test; end: # path(n, A) #-------------------------------------------------------------------------------- # Creates a path matrix of size (n + 1) by (n + 1), where 1's represent the # arrival/service sequence for a given row of the case matrix C. All other # elements in the path matrix are 0. The path starts at the lower left # corner of the matrix and moves to the upper right corner. The first leg of # the path is always the arrival of customer 1 represented by the entries in # the [n + 1, 1] and [n + 1, 2] positions. A 1 to the right of the previous # 1 signifies an arrival, while a 1 above the previous 1 signifies a service # completion. The procedure requires the arguments n, number of customers # and A, a row from the case matrix C. # # Name : path.mw # Author : Billy Kaczynski # Language : MAPLE 9 # Latest Revision : 01/12/09 # ---------------------------------------------------------------------------------- path := proc(n, A) local j, row, col, pat; row := n + 1; col := 2; pat := Matrix(n + 1, n + 1); pat[n + 1, 1] := 1; pat[n + 1, 2] := 1; for j from 2 to (2 * n) do if (A[1, j] = 1) then col := col + 1; pat[row, col] := 1; else row := row - 1; pat[row, col] := 1; fi; od: return pat; end: # kpath(n, k, A) #-------------------------------------------------------------------------------- # Creates a path matrix of size (n + k + 1) by (n + 1), where 1's represent the # arrival/service sequence for a given row of the case matrix C. All other # elements in the path matrix are 0. The path starts at the lower left # corner of the matrix and moves to the upper right corner. The first leg of # the path is either the arrival of a customer represented by the entry in # the [n + k + 1, 2] or a departure represented by the entry in the [n + k, 1] # positions. A 1 to the right of the previous 1 signifies an arrival, # while a 1 above the previous 1 signifies a service completion. The # procedure requires the arguments n, the number of customers processing through # the system arriving after time 0, k, the number of customers present at time 0, # and A, a row from the case matrix C. # # Name : kpath.mw # Author : Billy Kaczynski # Language : MAPLE 9 # Latest Revision : 03/07/09 # ---------------------------------------------------------------------------------- kpath := proc(n, k, A) local j, row, col, pat; row := n + k + 1; col := 1; pat := Matrix(n + k + 1, n + 1); pat[n + k + 1, 1] := 1; for j from 1 to (2 * n + k) do if (A[1, j] = 1) then col := col + 1; pat[row, col] := 1; else row := row - 1; pat[row, col] := 1; fi; od: return pat; end: # Cprime(n, C) #-------------------------------------------------------------------------------- # Produces the matrix defined as C', that is the distribution segment matrix # where each row represents the distribution segments for the case represented # by the corresponding row in the case matrix C. The elements of C' are # limited to a 0, 1, and 2, where 0 implies no sojourn time distribution # segment due to an emptying of the system, 1 implies a competing risk of an # arrival or completion of service and is distributed # exponential(lambda + mu), and a 2 implies a service completion # distribution leg which is distributed exponential(mu). The matrix C' # has the same number of rows as C, and 2n - 1 columns. Cprime(n, C) # calls path(n, A) and uses the path matrix to determine the appropriate # probability distribution function segments. The procedure requires the # arguments n, number of customers and C, the case matrix. # # Name : Cprime.mw # Author : Billy Kaczynski # Language : MAPLE 9 # Latest Revision : 01/12/09 # ---------------------------------------------------------------------------------- Cprime := proc(n, C) local prime, i, pat, dist, j, row, col; prime := Matrix(RowDimension(C), 2 * n - 1); for i from 1 to RowDimension(C) do row := n + 1; col := 2; pat := path(n, C[[i], 1 .. -1]); dist := Matrix(1, 2 * n - 1); for j from 1 to (2 * n - 1) do if (pat[row - 1, col] = 1) and (col < n + 1) then row := row - 1; dist[1, j] := 1; elif (pat[row - 1, col] = 1) and (col = n + 1) then row := row - 1; dist[1, j] := 2; elif (pat[row, col + 1] = 1) and (row + col > n + 2) then col := col + 1; dist[1, j] := 1; else col := col + 1; dist[1, j] := 0; fi; od; prime[[i], 1 .. -1] := dist; od; return prime; end: # kCprime(n, k, C) #-------------------------------------------------------------------------------- # Produces the matrix defined as C', that is the distribution segment matrix # where each row represents the distribution segments for the case represented # by the corresponding row in the case matrix C. The elements of C' are # limited to a 0, 1, and 2, where 0 implies no sojourn time distribution # segment due to an emptying of the system, 1 implies a competing risk of an # arrival or completion of service and is distributed # exponential(lambda + mu), and a 2 implies a service completion # distribution leg which is distributed exponential(mu). The matrix C' # has the same number of rows as C, and 2(n + 1) + k columns. Cprime(n, k, C) # calls kpath(n, k, A) and uses the path matrix to determine the appropriate # probability distribution function segments. The procedure requires the # arguments n, number of customers arriving after time 0, k, the number of # customers present at time 0, and C, the case matrix. # # Name : kCprime.mw # Author : Billy Kaczynski # Language : MAPLE 9 # Latest Revision : 03/07/09 # ---------------------------------------------------------------------------------- kCprime := proc(n, k, C) local prime, i, pat, dist, j, row, col; prime := Matrix(RowDimension(C), 2 * n + k ); for i from 1 to RowDimension(C) do row := n + k + 1; col := 1; pat := kpath(n, k, C[[i], 1 .. -1]); dist := Matrix(1, 2 * n + k ); for j from 1 to (2 * n + k ) do if (pat[row - 1, col] = 1) and (col < n + 1) then row := row - 1; dist[1, j] := 1; elif (pat[row - 1, col] = 1) and (col = n + 1) then row := row - 1; dist[1, j] := 2; elif (pat[row, col + 1] = 1) and (row + col > n + 2) then col := col + 1; dist[1, j] := 1; else col := col + 1; dist[1, j] := 0; fi; od; prime[[i], 1 .. -1] := dist; od; return prime; end: # caseprob(n, P) #-------------------------------------------------------------------------------- # Computes the probability associated with a given row of the case matrix C # as represented by the path created by path(n, A). Similar to how C' # identifies the appropriate distribution segments along the path, # caseprob(n, P) identifies the appropriate probability for each leg of the # path based on whether a competing risk occurs. Requires the arguments # n, number of customers and P, the path of a given case. Returns the # probability of the case passed to the procedure. # # Name : caseprob.mw # Author : Billy Kaczynski # Language : MAPLE 9 # Latest Revision : 01/12/09 # ---------------------------------------------------------------------------------- caseprob := proc(n, P) global X, Y; local p, j, row, col; p := 1; row := n + 1; col := 2; for j from 1 to (2 * n - 1) do if (P[row - 1, col] = 1) and (col < n + 1) then row := row - 1; p := p * 1 / Mean(Y) / (1 / Mean(X) + 1 / Mean(Y)); elif (P[row - 1, col] = 1) and (col = n + 1) then row := row - 1; elif (P[row, col + 1] = 1) and (row + col > n + 2) then col := col + 1; p := p * 1 / Mean(X) / (1 / Mean(X) + 1 / Mean(Y)); else col := col + 1; fi; od: return p; end: # kcaseprob(n, k, P) #-------------------------------------------------------------------------------- # Computes the probability associated with a given row of the case matrix C # as represented by the path created by kpath(n, k, A). Similar to how C' # identifies the appropriate distribution segments along the path, # caseprob(n, k, P) identifies the appropriate probability for each leg of the # path based on whether a competing risk occurs. Requires the arguments # n, number of customers arriving after time 0, k, the number of customers # present at time 0, and P, the path of a given case. Returns the # probability of the case passed to the procedure. # # Name : kcaseprob.mw # Author : Billy Kaczynski # Language : MAPLE 9 # Latest Revision : 03/07/09 # ---------------------------------------------------------------------------------- kcaseprob := proc(n, k, P) global X, Y; local p, j, row, col; p := 1; row := n + k + 1; col := 1; for j from 1 to (2 * n + k) do if (P[row - 1, col] = 1) and (col < n + 1) then row := row - 1; p := p * 1 / Mean(Y) / (1 / Mean(X) + 1 / Mean(Y)); elif (P[row - 1, col] = 1) and (col = n + 1) then row := row - 1; elif (P[row, col + 1] = 1) and (row + col > n + 2) then col := col + 1; p := p * 1 / Mean(X) / (1 / Mean(X) + 1 / Mean(Y)); else col := col + 1; fi; od: return p; end: # probvec(n) #-------------------------------------------------------------------------------- # Uses the procedure caseprob(n, P) successively to build a vector of # probabilities, one for each case of the C matrix. This vector has length # (2n)! / n! / (n = 1)!, which the the n-th Catalan number. Requires the # argument n, the number of customers. # # Name : probvec.mw # Author : Billy Kaczynski # Language : MAPLE 9 # Latest Revision : 01/12/09 # ---------------------------------------------------------------------------------- probvec := proc(n) local i, p, c; c := (2 * n)! / n! / (n + 1)!; p := Vector(c); for i from 1 to c do p[i] := caseprob(n, path(n, cases(n)[[i], 1 .. -1])); od: return p; end: # kprobvec(n, k) #-------------------------------------------------------------------------------- # Uses the procedure kcaseprob(n, k, P) successively to build a vector of # probabilities, one for each case of the C matrix. This vector has length # C. Requires the arguments n, the number of customers arriving after time 0 # and k, the number of customers present at time 0. # # Name : kprobvec.mw # Author : Billy Kaczynski # Language : MAPLE 9 # Latest Revision : 03/07/09 # ---------------------------------------------------------------------------------- kprobvec := proc(n, k) local i, p, C; C := kcases(n, k); p := Vector(RowDimension(C)); for i from 1 to RowDimension(C) do p[i] := kcaseprob(n, k, kpath(n, k, C[[i], 1 .. -1])); od: return p; end: # Tmat(a, b, A) #-------------------------------------------------------------------------------- # Creates the 2 by 2 matrix for determining whether selected customer # sojourn times are independent. Also provides information on the required # distribution segments for calculating the joint distribution between two # customers. Requires the arguments a, the first customer of interest in the # system, b, the second customer of interest in the system, and A, a single # row of the case matrix C, representing a given case. It uses this row of C # to identify the start and finish indices for customers a and b. If these # indices overlap, sojourn times are dependent, if they do not overlap # the sojourn times are independent. # # Name : Tmat.mw # Author : Billy Kaczynski # Language : MAPLE 9 # Latest Revision : 01/12/09 # ---------------------------------------------------------------------------------- Tmat := proc(a, b, A) local sta, fina, stb, finb, indexa, indexb, i, T; indexa := 0; indexb := 0; for i from 1 to ColumnDimension(A) do if A[1, i] = 1 then indexa := indexa + 1; if indexa = a then sta := i fi: if indexa = b then stb := i fi: elif A[1, i] = -1 then indexb := indexb - 1; if indexb = -(a) then fina := i fi: if indexb = -(b) then finb := i fi: fi: od: T := Matrix(2, 2, [[sta, fina], [stb, finb]]); return T; end: # kTmat(a, b, k, A) #-------------------------------------------------------------------------------- # Creates the 2 by 2 matrix for determining whether selected customer # sojourn times are independent and whether the customer index is less # than or equal to the number of customers present at time 0, k. # Also provides information on the required distribution segments for # calculating the joint distribution between two customers. Requires # the arguments a, the first customer of interest in the system, b, # the second customer of interest in the system, k, number of customers # present at time 0, and A, a single row of the case matrix C, # representing a given case. It uses this row of C to identify the # start and finish indices for customers a and b. If these # indices overlap, sojourn times are dependent, if they do not overlap # the sojourn times are independent. # # Name : kTmat.mw # Author : Billy Kaczynski # Language : MAPLE 9 # Latest Revision : 03/08/09 # ---------------------------------------------------------------------------------- kTmat := proc(a, b, k, A) local sta, fina, stb, finb, indexa, indexb, i, T; indexa := 0; indexb := 0; if a <= k then sta := 1; else for i from 1 to ColumnDimension(A) do if A[1, i] = 1 then indexa := indexa + 1; if indexa = a - k then sta := i + 1 fi: fi: od: fi: indexa := 0; if b <= k then stb := 1; else for i from 1 to ColumnDimension(A) do if A[1, i] = 1 then indexa := indexa + 1; if indexa = b - k then stb := i + 1 fi: fi: od: fi: for i from 1 to ColumnDimension(A) do if A[1, i] = -1 then indexb := indexb - 1; if indexb = -(a) then fina := i fi: if indexb = -(b) then finb := i fi: fi: od: T := Matrix(2, 2, [[sta, fina], [stb, finb]]); return T; end: # inde(a, b, T, A) #-------------------------------------------------------------------------------- # Calculates the case-specific joint cumulative distribution function # for customers a and b whose sojourn times are independent by # multiplying the CDF's of each customer. The individual customer CDF's # are calculated by determining the type and number of distribution legs # using the arguments a, the first customer of interest, b, the second # customer of interest, T, the resulting matrix from the call Tmat(a, b, A), # and A, the row of C' associated with the specific case. The CDF forms # for each case arise from appropriately defined random variables in APPL. # The procedure returns the joint cumulative distribution function in a # vector of length two, where both elements are identical in order to match # the piecewise result for customers with dependent sojourn times. # # Name : inde.mw # Author : Billy Kaczynski # Language : MAPLE 9 # Latest Revision : 01/12/09 # ---------------------------------------------------------------------------------- inde := proc(a, b, T, A) options remember; global X, Y; local i, dist1, dist2, jcdf, expa, expb; expa := 0; expb := 0; for i from 1 to (T[1, 2] - T[1, 1]) do if A[1, T[1, 1] + i - 1] = 1 then expa := expa + 1; elif A[1, T[1, 1] + i - 1] = 2 then expb := expb + 1; fi: od: if expb = 0 then dist1 := ErlangRV(1 / Mean(X) + 1 / Mean(Y), expa) elif expa = 0 then dist1 := ErlangRV(1 / Mean(Y), expb) else dist1 := [[unapply(conv(expa, expb), t)], [0, infinity], ["Continuous", "PDF"]]; fi: expa := 0; expb := 0; for i from 1 to (T[2, 2] - T[2, 1]) do if A[1, T[2, 1] + i - 1] = 1 then expa := expa + 1; elif A[1, T[2, 1] + i - 1] = 2 then expb := expb + 1; fi: od: if expb = 0 then dist2 := ErlangRV(1 / Mean(X) + 1 / Mean(Y), expa) elif expa = 0 then dist2 := ErlangRV(1 / Mean(Y), expb) else dist2 := [[unapply(conv(expa, expb), t)], [0, infinity], ["Continuous", "PDF"]]; fi: jcdf := apply(op(CDF(dist1)[1]), t[a]) * apply(op(CDF(dist2)[1]), t[b]); return Matrix([jcdf, jcdf]); end: # kinde(a, b, k, T, A) #-------------------------------------------------------------------------------- # Calculates the case-specific joint cumulative distribution function # for customers a and b whose sojourn times are independent by # multiplying the CDF's of each customer. The individual customer CDF's # are calculated by determining the type and number of distribution legs # using the arguments a, the first customer of interest, b, the second # customer of interest, T, the resulting matrix from the call kTmat(a, b, k, A), # and A, the row of C' associated with the specific case. The CDF forms # for each case arise from appropriately defined random variables in APPL. # The procedure returns the joint cumulative distribution function in a # vector of length two, where both elements are identical in order to match # the piecewise result for customers with dependent sojourn times. # # Name : kinde.mw # Author : Billy Kaczynski # Language : MAPLE 9 # Latest Revision : 03/15/09 # ---------------------------------------------------------------------------------- kinde := proc(a, b, k, T, A) options remember; global X, Y; local i, dist1, dist2, jcdf, expa, expb; expa := 0; expb := 0; for i from 0 to (T[1, 2] - T[1, 1]) do if A[1, T[1, 1] + i] = 1 then expa := expa + 1; elif A[1, T[1, 1] + i] = 2 then expb := expb + 1; fi: od: if expb = 0 then dist1 := ErlangRV(1 / Mean(X) + 1 / Mean(Y), expa) elif expa = 0 then dist1 := ErlangRV(1 / Mean(Y), expb) else dist1 := [[unapply(conv(expa, expb), w)], [0, infinity], ["Continuous", "PDF"]]; fi: expa := 0; expb := 0; for i from 0 to (T[2, 2] - T[2, 1]) do if A[1, T[2, 1] + i] = 1 then expa := expa + 1; elif A[1, T[2, 1] + i] = 2 then expb := expb + 1; fi: od: if expb = 0 then dist2 := ErlangRV(1 / Mean(X) + 1 / Mean(Y), expa) elif expa = 0 then dist2 := ErlangRV(1 / Mean(Y), expb) else dist2 := [[unapply(conv(expa, expb), w)], [0, infinity], ["Continuous", "PDF"]]; fi: jcdf := apply(op(CDF(dist1)[1]), t[a]) * apply(op(CDF(dist2)[1]), t[b]); return Matrix([jcdf, jcdf]); end: # dep(a, b, T, A) #-------------------------------------------------------------------------------- # Calculates the case-specific joint cumulative distribution function # for customers a and b whose sojourn times are dependent by # conditioning on the overlap distribution segment(s). The customer # sojourn time segments are divided up into their associated independent # and dependent (overlap) portions. This amounts to three segments, # customer a's independent portion defined as dist1, customer b's # independent portion defined as dist2, and the dependent overlap portion # defined as dist3. The joint cumulative distribution function has two # pieces, for the cases when t[a] < t[b] and t[b] < t[a]. # # Name : dep.mw # Author : Billy Kaczynski # Language : MAPLE 9 # Latest Revision : 01/12/09 # ---------------------------------------------------------------------------------- dep := proc(a, b, T, A) options remember; global X, Y; local i, expa, expb, dist1, dist2, dist3, jcdftop, jcdfbot; expa := 0; expb := 0; for i from 1 to (T[2, 1] - T[1, 1]) do if A[1, T[1, 1] + i - 1] = 1 then expa := expa + 1; elif A[1, T[1, 1] + i - 1] = 2 then expb := expb + 1; fi: od: if expb = 0 then dist1 := ErlangRV(1 / Mean(X) + 1 / Mean(Y), expa) elif expa = 0 then dist1 := ErlangRV(1 / Mean(Y), expb) else dist1 := [[unapply(conv(expa, expb), t)], [0, infinity], ["Continuous", "PDF"]]; fi: expa := 0; expb := 0; for i from 1 to (T[2, 2] - T[1, 2]) do if A[1, T[1, 2] + i - 1] = 1 then expa := expa + 1; elif A[1, T[1, 2] + i - 1] = 2 then expb := expb + 1; fi: od: if expb = 0 then dist2 := ErlangRV(1 / Mean(X) + 1 / Mean(Y), expa) elif expa = 0 then dist2 := ErlangRV(1 / Mean(Y), expb) else dist2 := [[unapply(conv(expa, expb), t)], [0, infinity], ["Continuous", "PDF"]]; fi: expa := 0; expb := 0; for i from 1 to (T[1, 2] - T[2, 1]) do if A[1, T[2, 1] + i - 1] = 1 then expa := expa + 1; elif A[1, T[2, 1] + i - 1] = 2 then expb := expb + 1; fi: od: if expb = 0 then dist3 := ErlangRV(1 / Mean(X) + 1 / Mean(Y), expa) elif expa = 0 then dist3 := ErlangRV(1 / Mean(Y), expb) else dist3 := [[unapply(conv(expa, expb), t)], [0, infinity], ["Continuous", "PDF"]]; fi: jcdftop := int(apply(op(CDF(dist1)[1]), t[a] - y) * apply(op(CDF(dist2)[1]), t[b] - y) * apply(op(PDF(dist3)[1]), y), y = 0 .. t[a]); jcdfbot := int(apply(op(CDF(dist1)[1]), t[a] - y) * apply(op(CDF(dist2)[1]), t[b] - y) * apply(op(PDF(dist3)[1]), y), y = 0 .. t[b]); return Matrix([jcdftop, jcdfbot]); end: # kdep(a, b, k, T, A) #-------------------------------------------------------------------------------- # Calculates the case-specific joint cumulative distribution function # for customers a and b whose sojourn times are dependent by # conditioning on the overlap distribution segment(s). The customer # sojourn time segments are divided up into their associated independent # and dependent (overlap) portions. This amounts to three segments, # customer a's independent portion defined as dist1, customer b's # independent portion defined as dist2, and the dependent overlap portion # defined as dist3. The joint cumulative distribution function has two # pieces, for the cases when t[a] < t[b] and t[b] < t[a]. When a and b are # both less than or equal to k, only two distribution segemtns arise, dist1 # and dist2, because both sojourn times start at time 0. Therefore the # sojourn time T[b] > T[a], and the resulting joint cumulative distribution # function has only a single piece. # # Name : kdep.mw # Author : Billy Kaczynski # Language : MAPLE 9 # Latest Revision : 03/15/09 # ---------------------------------------------------------------------------------- kdep := proc(a, b, k, T, A) options remember; global X, Y; local i, expa, expb, dist1, dist2, dist3, jcdftop, jcdfbot; expa := 0; expb := 0; if ((a <= k) and (b <= k)) then for i from 0 to (T[1, 2] - T[1, 1]) do if A[1, T[1, 1] + i] = 1 then expa := expa + 1; elif A[1, T[1, 1] + i] = 2 then expb := expb + 1; fi: od: #fi: if expb = 0 then dist1 := ErlangRV(1 / Mean(X) + 1 / Mean(Y), expa) elif expa = 0 then dist1 := ErlangRV(1 / Mean(Y), expb) else dist1 := [[unapply(conv(expa, expb), w)], [0, infinity], ["Continuous", "PDF"]]; fi: expa := 0; expb := 0; for i from 1 to (T[2, 2] - T[1, 2]) do if A[1, T[1, 2] + i ] = 1 then expa := expa + 1; elif A[1, T[1, 2] + i ] = 2 then expb := expb + 1; fi: od: if expb = 0 then dist2 := ErlangRV(1 / Mean(X) + 1 / Mean(Y), expa) elif expa = 0 then dist2 := ErlangRV(1 / Mean(Y), expb) else dist2 := [[unapply(conv(expa, expb), w)], [0, infinity], ["Continuous", "PDF"]]; fi: jcdftop := simplify(int(int(apply(op(PDF(dist2)[1]), y) * apply(op(PDF(dist1)[1]), x), y = 0 .. (t[b] - x)), x = 0 .. t[a])); return Matrix([jcdftop]); else for i from 1 to (T[2, 1] - T[1, 1]) do if A[1, T[1, 1] + i - 1] = 1 then expa := expa + 1; elif A[1, T[1, 1] + i - 1] = 2 then expb := expb + 1; fi: od: fi: if expb = 0 then dist1 := ErlangRV(1 / Mean(X) + 1 / Mean(Y), expa) elif expa = 0 then dist1 := ErlangRV(1 / Mean(Y), expb) else dist1 := [[unapply(conv(expa, expb), w)], [0, infinity], ["Continuous", "PDF"]]; fi: expa := 0; expb := 0; for i from 1 to (T[2, 2] - T[1, 2]) do if A[1, T[1, 2] + i] = 1 then expa := expa + 1; elif A[1, T[1, 2] + i] = 2 then expb := expb + 1; fi: od: if expb = 0 then dist2 := ErlangRV(1 / Mean(X) + 1 / Mean(Y), expa) elif expa = 0 then dist2 := ErlangRV(1 / Mean(Y), expb) else dist2 := [[unapply(conv(expa, expb), w)], [0, infinity], ["Continuous", "PDF"]]; fi: expa := 0; expb := 0; for i from 0 to (T[1, 2] - T[2, 1]) do if A[1, T[2, 1] + i] = 1 then expa := expa + 1; elif A[1, T[2, 1] + i] = 2 then expb := expb + 1; fi: od: if expb = 0 then dist3 := ErlangRV(1 / Mean(X) + 1 / Mean(Y), expa) elif expa = 0 then dist3 := ErlangRV(1 / Mean(Y), expb) else dist3 := [[unapply(conv(expa, expb), w)], [0, infinity], ["Continuous", "PDF"]]; fi: jcdftop := int(apply(op(CDF(dist1)[1]), t[a] - y) * apply(op(CDF(dist2)[1]), t[b] - y) * apply(op(PDF(dist3)[1]), y), y = 0 .. t[a]); jcdfbot := int(apply(op(CDF(dist1)[1]), t[a] - y) * apply(op(CDF(dist2)[1]), t[b] - y) * apply(op(PDF(dist3)[1]), y), y = 0 .. t[b]); return Matrix([jcdftop, jcdfbot]); end: # jpdf(a, b, n) #-------------------------------------------------------------------------------- # Creates the case joint cumulative distribution functions in a matrix with # dimension (2n)! / n! / (n + 1)! by 2. Calls the procedure Tmat(a, b, A) # and depending on the structure returned, calls the procedures # inde(a, b, T, A) or dep(a, b, T, A) to generate the appropriate case-wise # joint cumulative distribution function. Requires arguments a, the # first customer of interest, b, the second customer of interest, and n, # the number of customers. # # Name : jpdf.mw # Author : Billy Kaczynski # Language : MAPLE 9 # Latest Revision : 01/13/09 # ---------------------------------------------------------------------------------- jpdf := proc(a, b, n) local C, i, c, dist, T; C := cases(n); c := (2 * n)! / n! / (n + 1)!; dist := Matrix(c, 2); for i from 1 to c do T := Tmat(a, b, C[[i], 1 .. -1]); if T[1, 2] < T[2, 1] then dist[[i], 1 .. -1] := inde(a, b, T, Cprime(n, C[[i], 1 .. -1])); else dist[[i], 1 .. -1] := dep(a, b, T, Cprime(n, C[[i], 1 .. -1])); fi: od: return dist; end: # kjcdf(a, b, k, n) #-------------------------------------------------------------------------------- # Creates the case joint cumulative distribution functions in a matrix by # calling the procedure kTmat(a, b, k, A) and depending on the structure # returned, calls the procedures kinde(a, b, k, T, A) or kdep(a, b, k, T, A) # to generate the appropriate case-wise joint cumulative distribution function. # Requires arguments a, the first customer of interest, b, the second customer # of interest, and n, the number of customers arriving after time 0, and k, # the number of customers present at time 0. # # Name : kjcdf.mw # Author : Billy Kaczynski # Language : MAPLE 9 # Latest Revision : 03/15/09 # ---------------------------------------------------------------------------------- kjcdf := proc(a, b, k, n) local C, i, c, dist, T; C := kcases(n, k); #c := (2 * n)! / n! / (n + 1)!; dist := Matrix(RowDimension(C), 2); for i from 1 to RowDimension(C) do T := kTmat(a, b, k, C[[i], 1 .. -1]); if T[1, 2] < T[2, 1] then dist[[i], 1 .. -1] := kinde(a, b, k, T, kCprime(n, k, C[[i], 1 .. -1])); else dist[[i], 1 .. -1] := kdep(a, b, k, T, kCprime(n, k, C[[i], 1 .. -1])); fi: od: return dist; end: # conv(m, n) #-------------------------------------------------------------------------------- # Sums the appropriate distribution segments for the independent and dependent # portions of customer sojourn times bypassing the required calls to # Convolution(X, Y) in APPL by rewriting the integral as sums. Saves # significant CPU time by recognizing that these segments can all be written # as a sum of Erlang random variables. Requires the arguments m, the number # of expon(lambda + mu) segments and n, the number of expon(mu) segments. # # Name : conv.mw # Author : Billy Kaczynski # Language : MAPLE 9 # Latest Revision : 02/2/09 # ---------------------------------------------------------------------------------- conv := proc(m, n) options remember; global X, Y; local f1, f2, r, x, t, pdf; f1 := unapply((-1) ^ r * (m - 1 + x)! * (t) ^ (m - 1 + x - r) / (m - 1 + x - r)! / (1 / Mean(Y) - (1 / Mean(X) + 1 / Mean(Y))) ^ (r + 1), r, x); f2 := unapply((-1) ^ x * (n - 1)! / (n - 1 - x)! / x! * w ^ (n - 1 - x) * exp((1 / Mean(Y) - (1 / Mean(X) + 1 / Mean(Y))) * t), x); pdf := (1 / Mean(X) + 1 / Mean(Y)) ^ m * (1 / Mean(Y)) ^ n * exp(-(1 / Mean(Y)) * w) / (m - 1)! / (n - 1)! * add(f2(x) * add(f1(r, x), r = 0 .. m - 1 + x), x = 0 .. n - 1); return simplify(subs(t = w, pdf) - subs(t = 0, pdf)); end: # Cov(X, Y, a, b) #-------------------------------------------------------------------------------- # Mixes the results returned by probvec(n) and jpdf(a, b, n) to compute the # joint cumulative distribution function encompassing all cases. # Differentiates the results to produce the piecewise joint probability # distribution function. Calls Queue(n, k, s) to find the appropriate # expected values for the customers of interest, then uses the expected # values along with the expected value E(T[a]T[b]) found using the joint # probability distribution function to compute the covariance as # Cov(T[a], T[b]) = E(T[a]T[b]) - E(T[a])E(T[b]). Requires the arguments # X, the distribution of time between arrivals in the APPL list-of-list # format, Y, the service time distribution in the list-of-list format, # a, first customer of interest, and b, second customer of interest (a < b). # # Name : cov.mw # Author : Billy Kaczynski # Language : MAPLE 9 # Latest Revision : 01/13/09 # ---------------------------------------------------------------------------------- Cov := proc(X, Y, a, b) global lambda1, lambda2; local JPDFMAT, PVEC, JPDF, fta, ftb, Etatb, Eta, Etb, Cov; JPDFMAT := jpdf(a, b, b); PVEC := probvec(b); JPDF := Transpose(JPDFMAT) . PVEC; fta := simplify(diff(diff(JPDF[1], t[a]), t[b])); ftb := simplify(diff(diff(JPDF[2], t[a]), t[b])); Etatb := int(int(t[a] * t[b] * fta, t[a] = 0 .. t[b]), t[b] = 0 .. infinity) + int(int(t[a] * t[b] * ftb, t[b] = 0 .. t[a]), t[a] = 0 .. infinity); Eta := Mean(Queue(X, Y, a, 0, 1)); Etb := Mean(Queue(X, Y, b, 0, 1)); Cov := Etatb - Eta * Etb; return Cov; end: # kCov(X, Y, a, b, n, k) #-------------------------------------------------------------------------------- # Mixes the results returned by kprobvec(n, k) and kjpdf(a, b, n, k) to compute the # joint cumulative distribution function encompassing all cases. # Differentiates the results to produce the piecewise joint probability # distribution function. Calls Queue(n, k, s) to find the appropriate # expected values for the customers of interest, then uses the expected # values along with the expected value E(T[a]T[b]) found using the joint # probability distribution function to compute the covariance as # Cov(T[a], T[b]) = E(T[a]T[b]) - E(T[a])E(T[b]). Requires the arguments # X, the distribution of time between arrivals in the APPL list-of-list # format, Y, the service time distribution in the list-of-list format, # a, first customer of interest, and b, second customer of interest (a < b), # n, the number of customers arriving after time 0, and k, the number of # customers present at time 0. # # Name : kcov.mw # Author : Billy Kaczynski # Language : MAPLE 9 # Latest Revision : 03/17/09 # -------------------------------------------------------------------------------- kCov := proc(X, Y, a, b, n, k) local JPDFMAT, PVEC, JPDF, fta, ftb, Etatb, Eta, Etb, Cov, fa, fb, ftab; JPDFMAT := kjcdf(a, b, k, n); PVEC := kprobvec(n, k); JPDF := Transpose(JPDFMAT) . PVEC; if ((a <= k) and (b <= k)) then ftab := simplify(diff(diff(JPDF[1], t[a]), t[b])); Etatb := int(int(t[a]*t[b]*ftab, t[a] = 0..t[b]), t[b] = 0..infinity); Eta := int(t[a] * int(ftab, t[b] = t[a]..infinity), t[a] = 0..infinity); Etb := int(t[b] * int(ftab, t[a] = 0..t[b]), t[b] = 0..infinity); Cov := Etatb - Eta * Etb; return Cov; else fta := simplify(diff(diff(JPDF[1], t[a]), t[b])); ftb := simplify(diff(diff(JPDF[2], t[a]), t[b])); Etatb := int(int(t[a] * t[b] * fta, t[a] = 0 .. t[b]), t[b] = 0 .. infinity) + int(int(t[a] * t[b] * ftb, t[b] = 0 .. t[a]), t[a] = 0 .. infinity); fa := simplify(int(fta, t[b] = t[a]..infinity) + int(ftb, t[b] = 0..t[a])); Eta := int(t[a] * fa, t[a] = 0 .. infinity); fb := simplify(int(fta, t[a] = 0 .. t[b]) + int(ftb, t[a] = t[b] .. infinity)); Etb := int(t[b] * fb, t[b] = 0 .. infinity); Cov := Etatb - Eta * Etb; return Cov; fi: end: