# 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;