# copyright Barry L. Nelson and Michael R. Taaffe 2000
with(linalg):
with(plots):
PhPhInfNetMoments := proc(A, lambda, B, mu, K, Pr, nodelist, lower, upper)
#
# this procedure returns a plot of the mean and variance
# of number in at each node for the [Ph_t/Ph_t/Infinity]^K network as given
# in Nelson and Taaffe (2000)
#
# initial conditions for the network are empty and idle
#
# this version is designed to work with Maple 6/7/8
#
# last update 8/1/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
# nodelist = list of nodes whose performance is to be plotted
# lower = lower time limit for plot
# upper = upper time limit for plot
#
#
# Output:
# odeplot of mean and variance of number in queue for
# lower < t < upper
#
# [1] = mean plot
# [2] = variance plot
local mA, mB, N, NN, k, i, j, g1, a, b, c, d, mdeqs, v, off, meanplots,
varianceplots, r, M, V;
#
# create index mapping function
g1 := proc(i,j)
RETURN(N||i||N||j);
end;
# first check if there is only a single node
if K = 1 then
RETURN(PhPhInfMoments(A, lambda, B, mu, lower, upper));
else
# else it is a network, so use the network code
#
# obtain procedure to evaluate moment differential equations
mdeqs := PhPhInfNetNumeric(A, lambda, B, mu, K, Pr):
#
# determine number of phases in phase representations
# and offsets for nodes
mA := rowdim(A) -1;
v := vector(K);
off := vector(K+1);
v[1] := rowdim(B[1]) - 1;
off[1] := 0;
mB := v[1];
for r from 2 to K do
v[r] := rowdim(B[r]) - 1;
off[r] := off[r-1] + v[r-1];
od;
mB := off[K] + v[K];
off[K+1] := mB;
# create the N vector
N := matrix(mB, 1, [N||(1..mB)]);
N := N(t);
# create the NN matrix
NN := vector(mA);
for k from 1 to mA do
NN[k] := map(cat, matrix(mB,mB,g1), "A"||k)(t);
od;
# create the list of desired plots
meanplots := [t, sum(N[a,1], a=1..mB)];
varianceplots := [t, sum( sum( sum(NN[b][c,d], b=1..mA) -
N[c,1]*N[d,1], c=1..mB), d=1..mB) ];
for i in nodelist do
meanplots := meanplots, [ t, sum(N[a,1], a=off[i]+1..off[i+1]) ];
varianceplots := varianceplots,
[t, sum( sum( sum(NN[b][c,d], b=1..mA) -
N[c,1]*N[d,1], c=off[i]+1..off[i+1]), d=off[i]+1..off[i+1]) ];
od;
# create the mean and variance plot
M := odeplot(mdeqs, [ meanplots ], lower..upper, labels=[t, EN], color=black):
V := odeplot(mdeqs, [ varianceplots ], lower..upper, labels=[t, VarN], color=black):
RETURN(M, V);
fi;
end;