# copyright Barry L. Nelson and Michael R. Taaffe 1999
with(linalg):
with(plots):
PhPhInfSojourn := proc(A, lambda, B, mu, arrivalT)
#
# this procedure returns the mean, variance and a plot structure of
# the cdf of the sojourn time of a virtual arrival at time arrivalT
# to the Ph_t/Ph_t/Infinity queue as given in Nelson and
# Taaffe (1999)
#
# last update 6/30/99
#
# Inputs:
# A = Markov chain transition matrix for the arrival process
# B = Markov chain transition matrix for the service process
# lambda = arrival rate vector
# mu = service rate vector
# arrivalT = arrival time of virtual arrival
#
#
# Output:
# [1] mean of sojourn time
# [2] variance of sojourn time
# [3] plot structure for the cdf of sojourn time
#
local mA, mB, N, k, i, j, g1, a, c, d, b, x,
InitP, InitN, InitNA, InitNN, mdeq, cdf, first,
second, Tinf;
#
# create index mapping function
g1 := proc(i,j)
RETURN(N.i.N.j);
end;
#
# determine number of phases in phase representations
mA := rowdim(A) -1;
mB := rowdim(B) -1;
# extract beta
b := transpose(submatrix(B, mB+1..mB+1, 1..mB));
# create the initial conditions
InitP := vector(mA, 0);
InitN := vector(mB);
for i from 1 to mB do
InitN[i] := b[i,1];
od;
InitNA := matrix(mA, mB, 0);
InitNN := vector(mA);
for i from 1 to mA do
InitNN[i] := matrix(mB, mB, 0);
od;
# create the numerical procedure
mdeq := PhPhInfNumericInit(A, lambda, B, mu, arrivalT, InitP, InitN, InitNA, InitNN);
#
# create the N vector for plotting
N := matrix(mB, 1, [N.(1..mB)]);
N := N(t);
#
# compute first and second moment
cdf := subs(mdeq, 1 - sum(N[a,1], a=1..mB)):
first := subs(mdeq, M1(t));
second := subs(mdeq, M2(t));
Tinf := arrivalT;
while cdf(Tinf) < 0.99999 do
Tinf := Tinf + 0.1
od;
RETURN(first(Tinf), second(Tinf)- first(Tinf)^2,
odeplot(mdeq, [t, 1 - sum(N[a,1], a=1..mB)], arrivalT..arrivalT + Tinf));
end;