# copyright Barry L. Nelson and Michael R. Taaffe 2000 with(linalg): MapPhInfSymbolic := proc(A, vA, lambda, B, mu) # # this procedure returns a symbolic representation # of the first and second moment differential # equations for a Map_t/Ph_t/Infinity queue as given # in Nelson and Taaffe (2000c) # # initial conditions for the queue are empty and idle # # this version is designed to work with Maple V release 5 # # last update 11/1/00 # (created from PhPhInfSymbolic, last update 7/26/00) # # Inputs: # A = Markov chain transition matrix for the arrival process # vA = number of absorbing states in Map representation # lambda = arrival rate vector # B = Markov chain transition matrix for the service process # mu = service rate vector # # 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 # 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, A2rs, One, at, NNsum, NAsum, CTP; # # 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) -vA; 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(0); N := N(t); N0 := N0(0); # 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)(0); od; NA := matrix(mA,mB,g2)(t); NA0 := matrix(mA,mB,g2)(0); # create pieces of A and B A1 := submatrix(A, 1..mA, 1..mA); A2 := submatrix(A, 1..mA, mA+1..mA+vA); 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+vA, 1..mA)); b := transpose(submatrix(B, mB+1..mB+1, 1..mB)); One := matrix(vA, 1, 1); A2rs := multiply(A2, One); # 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(A2rs), 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] = 1; for i from 2 to mA+mB do deq := deq, dallF[i,1] = RHSF[i,1]; fcns := fcns, fcn[i,1]; initial := initial, fcn0[i,1] = 0; od; # # create cross-product moment differential equations # auxiliary procedures first... NAsum := proc() local NAtemp, NAbeta, jj; NAtemp := matrix(mB, mB, 0); for jj from 1 to mA do NAbeta := multiply( b, submatrix(NA,jj..jj, 1..mB)) + transpose(multiply( b, submatrix(NA,jj..jj, 1..mB))); NAtemp := evalm(NAtemp + scalarmul(NAbeta, lambda[jj]*( multiply(submatrix(a, i..i, 1..vA),submatrix(A2t, 1..vA, jj..jj))[1,1] ) )); od; RETURN(evalm(transpose(NAtemp))); end; NNsum := proc() local NNtemp, jj; NNtemp := matrix(mB, mB, 0); for jj from 1 to mA do NNtemp := evalm(NNtemp + scalarmul(NN[jj], lambda[jj]*( multiply(submatrix(a, i..i, 1..vA),submatrix(A2t, 1..vA, jj..jj))[1,1] ) )); od; RETURN(evalm(NNtemp)); end; CTP := proc() local CTtemp, jj; CTtemp := 0; for jj from 1 to mA do CTtemp := CTtemp + P[jj,1]*lambda[jj]*( multiply( submatrix(a, i..i, 1..vA),submatrix(A2t, 1..vA, jj..jj) )[1,1]); od; RETURN(evalm(D2(scalarmul(transpose(b), CTtemp), mB))); end; # # now create the terms 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]) + evalm(sum(NN[j]*(lambda[j]*A1[j,i]), j=1..mA) ) + NNsum() + NAsum() ) od; # # now correct the diagonal terms of the RHS CT := evalm(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) + CTP() ) 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] = 0; od; od; od; # # create first-partial moment differential equations RHSP := evalm( - multiply(lambdaD, NA) + multiply(a, A2t, lambdaD, NA) + transpose(multiply( b, transpose(multiply(a, A2t, C)) )) + 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] = 0; od; od; RETURN(deq, initial); end;