VacationQueueMoments[Lv_,Fv_,Pivs_,Bs_,Ls_,Fs_,Pisv_,n_,sigma_]:=Module[
{S,R,G,Psi,WorkingPrec,getRc,CRPSolve,DRPSolve,SigmaStar,SigmaStarDiff,Az,Ac,AS,M,β,c,stvphase,stsphase,sigAc,X,sigdiff,normc,MPower,NN,fmoms,tmp,tag},
WorkingPrec=120;
getRc[A0_,A1_,A2_,Prec_]:=Module[{dim,IA1,I,H,L,RR,T,Ro,U,IU},
dim=Dimensions[A1][[1]];
IA1=Inverse[-A1];
I=IdentityMatrix[dim];
H=A0.IA1;
L=A2.IA1;
RR=H;
T=L;
Ro=ConstantArray[1,{dim,dim}];
While [Norm[RR-Ro,∞]>Prec,
Ro=RR;
U=H.L+L.H;
IU=Inverse [I-U];
H=H.H.IU;
L=L.L.IU;
RR=RR+H.T;
T=L.T;
];
Return[RR];
];
CRPSolve[mbe_] :=Module[{mx,isize, p0},
	mx=mbe;
	isize=Dimensions[mx][[1]];
	p0=Table[0,{isize}];
	p0[[1]]=1 ;
	Do[mx[[i,1]]=1,{i, isize}];
	mx= Transpose[mx];
	Return[LinearSolve[mx, p0]];
];
DRPSolve[mbe_] :=CRPSolve[mbe-IdentityMatrix[Dimensions[mbe][[1]]]];
S=N[getRc[Bs,Ls,Fs,(10^(-WorkingPrec))],WorkingPrec];
R=N[getRc[Fs,Ls,Bs,(10^(-WorkingPrec))],WorkingPrec];
G = Inverse[-Ls-R.Bs].Bs;
Psi=Inverse[-S.Fs-Ls-R.Bs];
Az[z_] :=N[Lv+z Fv,WorkingPrec];
AS =N[KroneckerProduct[Transpose[Lv],IdentityMatrix[Dimensions[S][[1]]]] + KroneckerProduct[Transpose[Fv],S],WorkingPrec];
SigmaStar[A_]:=Quiet[NIntegrate[sigma[x] MatrixExp[A x],{x,0,∞},WorkingPrecision->WorkingPrec]];
SigmaStarDiff[nn_,c_]:=Quiet[NIntegrate[sigma[x](D[MatrixExp[Az[z] x],{z,nn}]/.z->c),{x,0,∞},WorkingPrecision->WorkingPrec]];
(*1. calculate the steady state distribution at the beginning of the cycles*)
M=KroneckerProduct[Transpose[Pisv],IdentityMatrix[Dimensions[Lv][[1]]]].SigmaStar[N[KroneckerProduct[IdentityMatrix[Dimensions[G][[1]]],Lv],WorkingPrec]+KroneckerProduct[Transpose[G],Fv]].ArrayReshape[Transpose[Pivs],{Dimensions[Pivs][[1]]Dimensions[Pivs][[2]],1}];
M=Transpose[ArrayReshape[M,{Dimensions[Lv][[1]],Dimensions[Lv][[2]]}]];
β = DRPSolve[M];
(* 2. calculate the normalization constant *)
c=N[1-(10^(-WorkingPrec /4)),WorkingPrec];
Ac = Az[c];
sigAc = SigmaStar[Ac];
X = (Transpose[ArrayReshape[Pivs.Psi,{Dimensions[Pivs.Psi][[1]]Dimensions[Pivs.Psi][[2]],1}]].SigmaStar[AS].KroneckerProduct[Transpose[{β}],IdentityMatrix[Dimensions[S][[1]]]])[[1]];
stvphase = β.(IdentityMatrix[Dimensions[Lv][[1]]]-sigAc).Inverse[-Ac];
stsphase = (β.sigAc.Pivs.Psi-X).(Inverse[IdentityMatrix[Dimensions[R][[1]]]-c R]+Inverse[c IdentityMatrix[Dimensions[S][[1]]]-S].S);
normc=Total[stvphase]+Total[stsphase];
(* Calculate factorial moments of the queue length *)
MPower[X_,x_]:=If[x==0,IdentityMatrix[Dimensions[X]],MatrixPower[X,x]];
NN=Dimensions[Lv][[1]] Dimensions[S][[1]];
fmoms = ConstantArray[0,{n+1}];
sigdiff={};
Do[
AppendTo[sigdiff, β.SigmaStarDiff[i,c]/Factorial[i]];
,{i,0,n}];
Do[
Do[
tmp=sigdiff[[d+1]].Pivs.Psi;
If[d==0,tmp-=X];
tag=Factorial[i]Total[(tmp.(MPower[R,i-d].MPower[Inverse[IdentityMatrix[Dimensions[R][[1]]]-c R],i-d+1]+(((-1))^(i-d))MPower[Inverse[c IdentityMatrix[Dimensions[S][[1]]]-S],i-d+1].S))]/normc;
(*Print[Precision[tag]];*)
fmoms[[i+1]]+=tag;
,{d,0,i}];
Do[
tmp =- sigdiff[[d+1]];
If[d==0,tmp+=β];
tag=Factorial[i]/Factorial[i-d] Total[tmp.(D[Inverse[-Lv-z Fv],{z,i-d}]/.z->c)]/normc;
(*Print[Precision[tag]];*)
fmoms[[i+1]]+=tag;
,{d,0,i}];
,{i,0,n}];
Return[fmoms];
];

Lv={{-5,1},{2,-3}};
Fv={{3,1},{0,1}};
Pivs={{8/10,0,2/10},{0,7/10,3/10}};
Ls={{-8,1,0},{0,-5,2},{1,3,-11}};
Bs={{0,0,0},{1,0,1},{4,1,2}};
Fs={{2,1,4},{0,1,0},{0,0,0}};
Pisv={{1,0},{1,0},{1/10,9/10}};
(* Parameters of the vacation period *)
α={7/10,2/10,1/10};
A={{-2,2,0},{0,-2,2},{0,0,-2}};
a=-Total[A,{2}];
sigmaPH[x_]:=α.MatrixExp[A x].a;
moms=VacationQueueMoments[Lv,Fv,Pivs,Bs,Ls,Fs,Pisv,2,sigmaPH]

VacMeans={1/1000,1/2,1,3/2,2,5/2,3,7/2,4,9/2,5,11/2,6,13/2,7,15/2,8,17/2,9,19/2,10};
Res={};
Do[
sigmaUni[x_]:=If[x<2mean,1/(2 mean),0];
sigmaExp[x_]:=1/mean (E^(-x/mean));
α={1,0,0};A={{-3/mean,3/mean,0},{0,-3/mean,3/mean},{0,0,-3/mean}};a=-Total[A,{2}];
sigmaErl[x_]:=α.MatrixExp[A x].a;
sigmaWei[x_]:=(1/2 2!/mean)(((x 2!/mean))^(1/2-1))(E^(-(((x 2!/mean))^(1/2))));
mUni=VacationQueueMoments2[Lv,Fv,Pivs,Bs,Ls,Fs,Pisv,1,sigmaUni][[2]];
mExp=VacationQueueMoments2[Lv,Fv,Pivs,Bs,Ls,Fs,Pisv,1,sigmaExp][[2]];
mErl=VacationQueueMoments2[Lv,Fv,Pivs,Bs,Ls,Fs,Pisv,1,sigmaErl][[2]];
mWei=VacationQueueMoments2[Lv,Fv,Pivs,Bs,Ls,Fs,Pisv,1,sigmaWei][[2]];
AppendTo[Res,{mUni,mExp,mErl,mWei}];
Print[Res];
,{mean,VacMeans}];