# copyright Barry L. Nelson and Michael R. Taaffe 2000
with(linalg):
PhPhInfNetSymbolic := proc(A, lambda, B, mu, K, Pr)
#
# this procedure returns a symbolic representation
# of the first and second moment differential
# equations for a K-node network of Ph_t/Ph_t/Infinity queues
# with a single arrival process
#
# initial conditions for the queue are empty and idle
#
# this version is designed to work with Maple 6/7/8
#
# 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
#
local BB, Bmu, i, j, g, h, v, r, offi, offj, mB;
# first check to make sure it is not a single-node problem
if K = 1 then
RETURN(PhPhInfSymbolic(A, lambda, B, mu));
# network problem; first extract key dimensions
else
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));
RETURN(PhPhInfSymbolic(A, lambda, BB, Bmu));
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;