# copyright Barry L. Nelson and Michael R. Taaffe 1999 with(linalg): PhPhInfNumericInit := proc(A, lambda, B, mu, InitT, InitP, InitN, InitNA, InitNN) # # this procedure returns a Maple listprocedure for # evaluating the first and second moment differential # equations for a Ph_t/Ph_t/Infinity queue as given # in Section 4 of Nelson and Taaffe (1999) # # it also adds differential equations for the first and # second moment of virtual sojourn time which only make # sense when there is initially a single customer and # no new arrivals # # initial conditions for the queue are arbitrary # # this version is designed to work with Maple 6/7/8 # # last update 6/30/99 7/27/00 08/20/01 10/15/03 # # Inputs: # A = Markov chain transition matrix for the arrival process # lambda = arrival rate vector # B = Markov chain transition matrix for the service process # mu = service rate vector # InitT = time at which the initial conditions apply # InitP = vector of length mA with initial conditions for P # InitN = vector of length mB with initial conditions for N # InitNA = mA x mB matrix of initial conditions for NA # InitNN = mA matrices of mB x mB with initial conditions for NN # # Outputs: # Pk = probability of being in arrival phase k # Nj = expected number in service phase j # NiNjAk = partial expected product of NiNj with respect to service phase is k # NjAk = partial expected value of Nj with respect to service phase k # M1 = first moment of sojourn time # M2 = second moment of sojourn time local A1, B1, A2, B2, a, b, C, D, mA, mB, IA, IB, Ep, Pp, RHS, all, deq, i, dall, N, P, NN, NA, muD, lambdaD, A2t, j, k, g1, g2, D2, D1, Diag, N0, P0, RHSC, RHSF, RHSP, dallC, dallF, dallP, fcn, fcn0, fcns, initial, NN0, NA0, CT; # # create index mapping functions g1 := proc(i,j) RETURN(N||i||N||j); end; g2 := proc(k,j) RETURN(N||j||"A"||k); end; D2 := proc(x, d) local i, M; M := matrix(d,d,0); for i from 1 to d do M[i,i] := x[1,i] od; RETURN(evalm(M)); end; D1 := proc(x, d) local i, M; M := matrix(d,d,0); for i from 1 to d do M[i,i] := x[i] od; RETURN(evalm(M)); end; Diag := proc(x) local i, k, M; k := rowdim(x); M := matrix(k,k,0); for i from 1 to k do M[i,i] := x[i,i] od; RETURN(M); end; # # determine number of phases in phase representations mA := rowdim(A) -1; mB := rowdim(B) -1; # create the N and P vectors P := matrix(mA, 1, [P||(1..mA)]); N := matrix(mB, 1, [N||(1..mB)]); P0 := matrix(mA, 1, [P||(1..mA)]); N0 := matrix(mB, 1, [N||(1..mB)]); P := P(t); P0 := P0(InitT); N := N(t); N0 := N0(InitT); # create the NN and NA matrices NN := vector(mA); NN0 := vector(mA); for k from 1 to mA do NN[k] := map(cat, matrix(mB,mB,g1), "A"||k)(t); NN0[k] := map(cat, matrix(mB,mB,g1), "A"||k)(InitT); od; NA := matrix(mA,mB,g2)(t); NA0 := matrix(mA,mB,g2)(InitT); # create pieces of A and B A1 := submatrix(A, 1..mA, 1..mA); A2 := submatrix(A, 1..mA, mA+1..mA+1); A2t := transpose(A2); B1 := submatrix(B, 1..mB, 1..mB); B2 := submatrix(B, 1..mB, mB+1..mB+1); a := transpose(submatrix(A, mA+1..mA+1, 1..mA)); b := transpose(submatrix(B, mB+1..mB+1, 1..mB)); # create diagonal rate matrices muD := D1(mu, mB); lambdaD := D1(lambda, mA); # create auxiliary matrices C and D C := multiply(lambdaD, P); D := multiply(muD, N); # form the first-moment differential equations IA := band([1], mA); IB := band([1], mB); Ep := multiply(b, multiply(transpose(A2), C)) + multiply(transpose(B1) - IB, D); Pp := multiply(a, multiply(transpose(A2), C)) + multiply(transpose(A1) - IA, C); RHSF := stackmatrix(Pp, Ep); dallF := map(diff, stackmatrix(P, N), t); fcn := stackmatrix(P, N); fcn0 := stackmatrix(P0, N0); deq := dallF[1,1] = RHSF[1,1]; fcns := fcn[1,1]; initial := fcn0[1,1] = InitP[1]; for i from 2 to mA do deq := deq, dallF[i,1] = RHSF[i,1]; fcns := fcns, fcn[i,1]; initial := initial, fcn0[i,1] = InitP[i]; od; for i from mA+1 to mA + mB do deq := deq, dallF[i,1] = RHSF[i,1]; fcns := fcns, fcn[i,1]; initial := initial, fcn0[i,1] = InitN[i-mA]; od; # # create cross-product moment differential equations RHSC := vector(mA); for i from 1 to mA do RHSC[i] := evalm(multiply(transpose(B1), muD, NN[i]) + multiply(NN[i], muD, B1) - multiply(muD, D2(submatrix(NA,i..i,1..mB),mB), B1) - multiply(transpose(B1), muD, D2(submatrix(NA,i..i,1..mB),mB)) - multiply(muD, NN[i]) - multiply(NN[i], muD) - scalarmul(NN[i], lambda[i]) + scalarmul(multiply(b, transpose(A2), lambdaD, NA), a[i,1]) + scalarmul( multiply(transpose(NA), lambdaD, A2, transpose(b)),a[i,1]) + evalm(sum(NN[j]*(lambda[j]*(A1[j,i] + a[i,1]*A2[j,1])), j=1..mA))) od; # # now correct the diagonal terms of the RHS CT := evalm( scalarmul(multiply(a, transpose(b)), multiply(A2t,C)[1,1]) + multiply(NA, muD) + multiply(NA, muD, B1) ); for i from 1 to mA do RHSC[i] := evalm( RHSC[i] + D2(submatrix(CT, i..i, 1..mB), mB) ) od; # # finish cross-product differential equations dallC := vector(mA); for i from 1 to mA do dallC[i] := map(diff, NN[i], t) od; for i from 1 to mA do for j from 1 to mB do for k from 1 to mB do deq := deq, dallC[i][j,k] = RHSC[i][j,k]; fcns := fcns, NN[i][j,k]; initial := initial, NN0[i][j,k] = InitNN[i][j,k]; od; od; od; # # create first-partial moment differential equations RHSP := evalm( - multiply(lambdaD, NA) + multiply(a, A2t, lambdaD, NA) + evalm(multiply(a, transpose(b))*multiply(A2t, C)[1,1]) + multiply(transpose(A1), lambdaD, NA) + multiply(NA, muD, B1) - multiply(NA, muD)); dallP := map(diff, NA, t); for i from 1 to mA do for j from 1 to mB do deq := deq, dallP[i,j] = RHSP[i,j]; fcns := fcns, NA[i,j]; initial := initial, NA0[i,j] = InitNA[i,j]; od; od; # # add differential equation for 1st & 2nd moment of sojourn time deq := deq, diff(M1(t), t) = sum(N[z,1], z=1..mB); fcns := fcns, M1(t); initial := initial, M1(InitT) = 0; deq := deq, diff(M2(t), t) = 2*(t - InitT)*sum(N[z,1], z=1..mB); fcns := fcns, M2(t); initial := initial, M2(InitT) = 0; RETURN(dsolve({deq, initial}, {fcns}, type=numeric, output=listprocedure)); end;