# copyright Barry L. Nelson and Michael R. Taaffe 2000
with(linalg):
with(plots):
PhPhInfNetSojourn := proc(A, lambda, B, mu, K, Pr, 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]^K network as given in Nelson and
# Taaffe (2000)
#
# last update 7/28/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
# [vector of matrices of dimension K]
# mu = service rate vector
# [vector of vectors of dimension K]
# K = number of nodes in the network
# Pr = Markov routing matrix
# 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, BB, Bmu, offi, offj, v, g, h, r;
#
# create index mapping function
g1 := proc(i,j)
RETURN(N||i||N||j);
end;
# check if there is only a single node
if K = 1 then
Return(PhPhInfSojourn(A, lambda, B, mu, arrivalT));
else
# else it is a network, so use the network code
#
# determine number of phases in phase representations
# and offsets for nodes
mA := rowdim(A) -1;
# first extract key dimensions
v := vector(K);
v[1] := rowdim(B[1]) - 1;
mB := v[1];
for r from 2 to K do
v[r] := rowdim(B[r]) - 1;
mB := mB + v[r];
od;
BB := matrix(mB+1, mB+1, 0);
# build the composite service-rate vector
Bmu := vector(mB);
offi := 0;
for i from 1 to K do
for g from 1 to v[i] do
Bmu[offi + g] := mu[i][g];
od;
offi := offi + v[i]
od;
# build the composite phase transition matrix
offi := 0; offj := 0;
# build rows and columns 1 through mB
for i from 1 to K do
for j from 1 to K do
for g from 1 to v[i] do
for h from 1 to v[j] do
BB[offi + g, offj + h] := Pr[i,j]*B[i][g, v[i]+1]*B[j][v[j]+1,h];
od;
od;
offj := offj + v[j];
od;
offi := offi + v[i];
offj := 0;
od;
# build the last row
offj := 0;
for j from 1 to K do
for h from 1 to v[j] do
BB[mB+1, offj + h] := Pr[K+1,j]*B[j][v[j]+1,h];
od;
offj := offj + v[j];
od;
# build the last column
offi := 0;
for i from 1 to K do
for g from 1 to v[i] do
BB[offi + g, mB+1] := Pr[i, K+1]*B[i][g, v[i]+1];
od;
offi := offi + v[i];
od;
BB[mB+1, mB+1] := 0;
# add the terms representing immediate feedback
BB := matadd(BB, MatDiag(B, K));
# extract beta
b := transpose(submatrix(BB, 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, BB, Bmu, 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));
fi;
end;
MatDiag := proc(x, k)
# auxiliary proc to build a matrix with other matrices on the diagonal
local m, lower, i, j, l, d, M;
m := rowdim(x[1]) - 1;
for l from 2 to k do
m := m + rowdim(x[l]) - 1;
od;
M := matrix(m+1, m+1, 0);
lower := 0;
for l from 1 to k do
d := rowdim(x[l]) - 1;
for i from 1 to d do
for j from 1 to d do
M[lower + i, lower + j] := x[l][i,j];
od;
od;
lower := lower + d;
od;
RETURN(M);
end;