# copyright Barry L. Nelson and Michael R. Taaffe 1999

with(linalg):

PhPhInfNumericEmpty  := proc(A, lambda, B, mu)
#
# 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) 
# 
# initial conditions for the queue are empty and idle
#
# 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
#
# 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
#

           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(0);
          N := N(t);
          N0 := N0(0);
# 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)(0);
          od;
          NA := matrix(mA,mB,g2)(t);
          NA0 := matrix(mA,mB,g2)(0);
# 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] = 1;

          for i from 2 to mA+mB do
             deq := deq, dallF[i,1] = RHSF[i,1];
             fcns := fcns, fcn[i,1];
             initial := initial, fcn0[i,1] = 0;
          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] = 0;
           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] = 0;
       od;
   od;
RETURN(dsolve({deq, initial}, {fcns}, type=numeric, method=rkf45, abserr=.0001, maxfun=-1,
            relerr = .1e-3, output=listprocedure));
end;