待ち行列の理論 本間鶴千代先生 理工学社

本間鶴千代著「待ち行列の理論」をプログラムしてみました。

# 待ち行列の理論

# M/M/s　p.21

library(dplyr)
library(pforeach)

# customer arrived data

total_customer_num=200
window_num=2
date=Sys.time()

cus_data=data.frame(customer_id=1:total_customer_num,arrived_time=0,cus_num=1) 

cus_data$arrived_time[1]=rexp(1,rate=0.07);cus_data$cus_num[1]=1
for(j in 2:nrow(cus_data)){
  
  cus_data$arrived_time[j]=cus_data$arrived_time[j-1]+sum(rexp(1,rate=0.07))
}

arrived_data=cus_data 

arrived_time=cus_data %>% dplyr::select(customer_id,arrived_time)

arrived_time=arrived_time[order(arrived_time$arrived_time),]

arrived_time=arrived_time %>% dplyr::mutate(diff_time=0)

for(j in 2:nrow(arrived_time)){
  arrived_time$diff_time[j]=arrived_time$arrived_time[j]-arrived_time$arrived_time[j-1]
}



# service time data 

cus_data=cus_data %>% dplyr::mutate(service_time=0)

for(j in 1:nrow(cus_data)){
  cus_data$service_time[j]=sum(rexp(cus_data$cus_num[j],rate = 0.08))
}

service_time=cus_data %>% dplyr::select(customer_id,cus_num,service_time) 

# analysis for M/M/S


# 到着時間のデータを指数分布で最尤推定量で当てはめる

lambda=1/mean(arrived_time$diff_time)

time=rexp(nrow(arrived_time),rate=lambda)

arrived_time=arrived_time %>% dplyr::mutate(predict=time)

plot(seq(1,nrow(arrived_time),1),arrived_time$diff_time,type="p",xlab="標本数",main="到着した後に次にお客様が到着するまでの時間",col=2,ylab="時間の長さ",ylim=c(min(time,arrived_time$diff_time),max(time,arrived_time$diff_time)))
par(new=T)
plot(seq(1,nrow(arrived_time),1),time,type="p",xlab="標本数",main="到着した後に次にお客様が到着するまでの時間",col=3,ylab="時間の長さ",ylim=c(min(time,arrived_time$diff_time),max(time,arrived_time$diff_time)))

print(paste0("サンプルの平均=",mean(arrived_time$diff_time),";サンプルの分散=",var(arrived_time$diff_time)))
print(paste0("予測値の平均=",mean(time),";予測値の分散=",var(time)))

# サービス時間のデータを指数分布で最尤推定量で当てはめる

mu=1/mean(service_time$service_time)

time=rexp(nrow(service_time),rate=mu)

service_time=service_time %>% dplyr::mutate(predict=time)

plot(seq(1,nrow(service_time),1),service_time$service_time,xlab="標本数",ylab="時間の長さ",main="お客様一人当たりに対するサービス時間",type="p",col=2,ylim=c(min(time,service_time$service_time),max(time,service_time$service_time)))
par(new=T)
plot(seq(1,nrow(service_time),1),time,xlab="標本数",ylab="時間の長さ",main="お客様一人当たりに対するサービス時間",type="p",col=3,ylim=c(min(time,service_time$service_time),max(time,service_time$service_time)))
print(paste0("サンプルの平均=",mean(service_time$service_time),";サンプルの分散=",var(service_time$service_time)))
print(paste0("予測値の平均=",mean(time),";予測値の分散=",var(time)))


# 平衡状態の選定

p=lambda/(mu*window_num)

if(p<1){
  print(paste0("p値は",p,"で平衡状態である。"))
  }else{
    print("平衡状態でない。")
  }
  
# 平衡状態における列の長さの確率

a=lambda/mu;s=window_num;

row_num_data=data.frame(row_num=0:1000,prob=0)

val=0
for(j in 0:(s-1)){
  val=val+(a^j)/(gamma(j+1))
}
row_num_data$prob[1]=p0=1/(sum(val+(a^s)/(gamma(s)*(s-a))))

for(j in 2:nrow(row_num_data)){
  if(j < s){
  row_num_data$prob[j]=p0*(a^j)/(gamma(j+1))
  }else{
  row_num_data$prob[j]=((s^s)*(p^j)*p0)/gamma(s+1)
  }
}

# 窓口がすべてふさがっている確率

prob=1-sum(row_num_data$prob[row_num_data$row_num<s])

print(paste0("窓口がすべてふさがっている確率は",prob))

# 列の長さ

L=((s^s)*(p^(s+1))*(p0))/(gamma(s+1)*(1-p)^2)

print(paste0("列の長さ平均は",L))

# 列の長さとサービス中人数の和

total_L=L+a-a*(1-sum(row_num_data$prob[row_num_data$row_num<s]))+s*(1-sum(row_num_data$prob[row_num_data$row_num<(s+1)]))

print(paste0("列の長さとサービス中の人数の平均は",total_L))

# waiting timeにおける平均待ち時間

waiting_time_data=data.frame(waiting_customers=1:10,ave_wait_time=0)

for(j in 1:nrow(waiting_time_data)){
  waiting_time_data$ave_wait_time[j]=(j+1)/(s*mu)
}

# 待ち時間の平均

W_q=((s^(s-1))*(p^s)*p0)/(gamma(s+1)*mu*(1-p)^2)

# W_q時間おきに来る平均訪問来客数

ave_cus=lambda*W_q

if(lambda*W_q==L){
  print("result is right!")
}

# t時間より多く待たされる確率 P(T>t)

wait_prob_data=data.frame(t=0:20,prob=0)
for(j in 1:nrow(wait_prob_data)){
wait_prob_data$prob[j]=(exp(-(1-p)*s*mu*wait_prob_data$t[j])*row_num_data$prob[row_num_data$row_num==s])/(1-p)

}

plot(wait_prob_data$t,wait_prob_data$prob,type="p",col=2,xlab="time",ylab="prob",main="待ち時刻の長さによって発生する確率")
  
print(paste0("待たされる確率は",wait_prob_data$prob[1]))

# G/M/1 p.38

# a(t)が指数分布の場合(M/M/1)；

library(dplyr)

# customer arrived data

total_customer_num=200

window_num=1 #固定

date=Sys.time()

cus_data=data.frame(customer_id=1:total_customer_num,arrived_time=0,cus_num=0) 

cus_data$arrived_time[1]=rexp(1,rate=0.04);cus_data$cus_num[1]=1

for(j in 2:nrow(cus_data)){
  
  cus_data$arrived_time[j]=cus_data$arrived_time[j-1]+sum(rexp(1,rate=0.04))
    
  cus_data$cus_num[j]=1

  }

arrived_time=cus_data %>% dplyr::select(customer_id,arrived_time)

arrived_time=arrived_time[order(arrived_time$arrived_time),]

arrived_time=arrived_time %>% dplyr::mutate(diff_time=0)

for(j in 2:nrow(arrived_time)){

    arrived_time$diff_time[j]=arrived_time$arrived_time[j]-arrived_time$arrived_time[j-1]

    }

# service time data 

cus_data=cus_data %>% dplyr::mutate(service_time=0)

for(j in 1:nrow(cus_data)){

 cus_data$service_time[j]=rexp(1,rate = 0.05)

}

service_time=cus_data %>% dplyr::select(customer_id,cus_num,service_time) %>% dplyr::mutate(customer_id,service_time=service_time)

# analysis for M/M/1

# 到着時間のデータを指数分布で最尤推定量で当てはめる

lambda=1/mean(arrived_time$diff_time)

time=rexp(nrow(arrived_time),rate=lambda)

arrived_time=arrived_time %>% dplyr::mutate(predict=time)

plot(seq(1,nrow(arrived_time),1),arrived_time$diff_time,type="p",xlab="標本数",main="到着した後に次にお客様が到着するまでの時間",col=2,ylab="時間の長さ",ylim=c(min(time,arrived_time$diff_time),max(time,arrived_time$diff_time)))

par(new=T)

plot(seq(1,nrow(arrived_time),1),time,type="p",xlab="標本数",main="到着した後に次にお客様が到着するまでの時間",col=3,ylab="時間の長さ",ylim=c(min(time,arrived_time$diff_time),max(time,arrived_time$diff_time)))

print(paste0("サンプルの平均=",mean(arrived_time$diff_time),";サンプルの分散=",var(arrived_time$diff_time)))

print(paste0("予測値の平均=",mean(time),";予測値の分散=",var(time)))

# サービス時間のデータを指数分布で最尤推定量で当てはめる

mu=1/mean(service_time$service_time)

time=rexp(nrow(service_time),rate=mu)

service_time=service_time %>% dplyr::mutate(predict=time)

plot(seq(1,nrow(service_time),1),service_time$service_time,xlab="標本数",ylab="時間の長さ",main="お客様一人当たりに対するサービス時間",type="p",col=2,ylim=c(min(time,service_time$service_time),max(time,service_time$service_time)))

par(new=T)

plot(seq(1,nrow(service_time),1),time,xlab="標本数",ylab="時間の長さ",main="お客様一人当たりに対するサービス時間",type="p",col=3,ylim=c(min(time,service_time$service_time),max(time,service_time$service_time)))

print(paste0("サンプルの平均=",mean(service_time$service_time),";サンプルの分散=",var(service_time$service_time)))

print(paste0("予測値の平均=",mean(time),";予測値の分散=",var(time)))


# 平衡状態の選定

p=lambda/(mu*window_num)

if(p<1){
 
print(paste0("p値は",p,"で平衡状態である。"))
  
  }else{
  
print("平衡状態でない。")
  
    }

z=1;y=1;v=(1/(lambda*mu))*(lambda+mu+z);

# y*x^2+lambda*v*x+p=0の解を求めることによって、モーメントを計算する

r=sqrt(v^2-(4*y)/(lambda*mu))

ans1=(lambda*(v-r))/(2*y)

ans2=(lambda*(v+r))/(2*y)

if(abs(ans1)<1){

  ans=ans1

}else{

  ans=ans2

}

p_pi=function(y,z){

w=(y*(1-ans))/(z/mu+1-y*ans)

return(w)

}

# 微分して係数をy=1,z=0とすることによってNを計算する。Nはサービス時間中にサービスを受け終わる平均人数となる。

h=0.0001

N=(p_pi(1+h,0)-p_pi(1-h,0))/(2*h)

for(m in 1:56){


Prob=function(t){
  
  if(m==1){
    
  w=((1/lambda)*(exp(-(lambda+mu)*t))*t^(2*m-2))/((1/((lambda*mu)^m))*gamma(m)*gamma(m+1))
  
  }else{
  
  w=((1/lambda)*(exp(-(lambda+mu)*t*(1/(2*m-2)))*t)^(2*m-2))/((1/((lambda*mu)^m))*gamma(m)*gamma(m+1))
  
  }

  return(w)

}

png(paste0("~/待ち行列(サービスを受け終わる人数の確率密度)/ サービスを受け終わる人数が",m,"である確率密度.png"))


plot(seq(0,5000,0.1),Prob(seq(0,5000,0.1)),type="p",col=2,xlab="時刻",ylab="確率",main=paste0("サービスを受け終わる人数が",m,"である確率密度"))

dev.off()

plot(seq(0,5000,0.1),Prob(seq(0,5000,0.1)),type="p",col=2,xlab="時刻",ylab="確率",main=paste0("サービスを受け終わる人数が",m,"である確率密度"))

}

# 修理工数の決定方法

# カメラの修理屋さんの例で、1時間当たりの修理要求台数のデータを作成

camera_repair=data.frame(num=2:16,freq=c(1,3,5,8,12,14,14,13,10,8,5,3,2,1,1))

# ポアソン分布で最尤推定

lambda=sum(camera_repair$num*(camera_repair$freq/100))

# 修理時間(単位時間1時間として0.25時間等を考える)

repair_time=data.frame(time=0.25*c(1:14),freq=c(26,19,14,10,7,7,6,3,2,2,1,1,1,1))

mu=1/sum(repair_time$time*(repair_time$freq/100))

a=lambda/mu



print(paste0("修理待ちのカメラをため込まないためには、少なくとも",ifelse(a-floor(a)>0,floor(a)+1,floor(a)),"人必要"))

library(dplyr)

# 機械修理の問題　p.58

machine_num=m=20

# 修理項数

repair_man=s=11

n=1000

# 修理時刻のデータ

rest_time=rexp(n,rate=0.06)

# 指数分布でモデルをフィットさせる
mu=1/mean(rest_time)

# 平均稼働時間のデータ

working_time_data=rexp(n,rate=0.05)

# 指数分布でモデルを合わせる

lambda=1/mean(working_time_data)

# 平衡状態の確認

p=lambda/(mu*s)

if(p<1){
 
   print(paste0("p値は",p,"で平衡状態である。"))
  
  }else{
  
      print("平衡状態でない。")
  }

# 系の長さの確率(修理を待っている機械の数の確率)

prob=data.frame(num=0:m,prob=0)

a=(lambda/mu)

p0=0

p0_1=0; p0_2=0;

for(j in 0:s){
  
  if(j==0){
 
  p0_1=p0_1+1*(a^j) 

  }else{

  p0_1=p0_1+prod(c(1:m)/c(c(1:j),c(1:(m-j))))*(a^j)


    }

  }

for(j in (s+1):m){
  
  if(j==m){
 
    p0_2=p0_2+1*(a^j)*(1/(prod(c(1:s))*(s^(j-s))))*prod(c(1:j))
  
  }else{
  
  p0_2=p0_2+prod(c(1:m)/c(c(1:j),c(1:(m-j))))*prod(c(1:j))*(a^j)*(1/(prod(c(1:s))*(s^(j-s))))
  
}

}

p0=1/(p0_1+p0_2)

for(j in 1:nrow(prob)){
 
  k=prob$num[j]
  
  if(s>k){
  
  if(k==0){
  
  prob$prob[j]=(a^k)*p0 
  
   }else{
  
  prob$prob[j]=prod(c(1:m)/c(c(1:k),c(1:(m-k))))*(a^k)*p0
  
  }
  
  }else{
  
  if(k==m){
  
    prob$prob[j]=((a^k)*p0*prod(c(1:k)))/(prod(c(1:s))*s^(k-s))
  
  }else{
  
  prob$prob[j]=(prod(c(1:m)/c(c(1:k),c(1:(m-k))))*(a^k)*p0*prod(c(1:k)))/(prod(c(1:s))*s^(k-s))
  
}

  }

}

plot(prob$num,prob$prob,xlab="修理を待っている機械の数",ylab="確率",main=paste0("修理中の機械の数の確率(修理工の数が",s,"の時)"),type="o",col=2)

# 修理を待っている機械の平均個数

L=sum((prob$num[prob$num>s]-s)*prob$prob[prob$num>s])

print(paste0("修理を待っている機械の平均個数は",L))

print(paste0("稼働機械数の平均は",sum((m-prob$num)*prob$prob)))

print(paste0("修理中の機械数の平均は",sum((prob$num[prob$num<s])*prob$prob[prob$num<s])+sum(s*prob$prob[prob$num>=s])))

library(dplyr)

# 予備品数の決定問題 p.62

# 公衆電話機が1日平均1台故障し、修理に1週間かかるとする。この時、予備として何台仕入れておけばよいか。

lambda=3 #1日。指数分布に従うとする。

mu=1/7 #1週間に一度損傷が起こると仮定

a=lambda/mu

pai=function(s){
  
p=lambda/(mu*s)

# 平衡状態における列の長さの確率

row_num_data=data.frame(row_num=0:1000,prob=0)

val=0

for(j in 0:(s-1)){

    val=val+(a^j)/(gamma(j+1))
}

row_num_data$prob[1]=p0=1/(sum(val+(a^s)/(gamma(s)*(s-a))))

for(j in 2:nrow(row_num_data)){

    if(j < s){
  
      row_num_data$prob[j]=p0*(a^j)/(gamma(j+1))
  
      }else{
  
        row_num_data$prob[j]=((s^s)*(p^j)*p0)/gamma(s+1)
  
        }

  }

y<-1-sum(row_num_data$prob[row_num_data$row_num<s])

return(y)

}

prob_data=data.frame(num=(round(a)+1):30,prob=0)

for(j in 1:nrow(prob_data)){

    prob_data$prob[j]=pai(prob_data$num[j])

    }

plot(prob_data$num,prob_data$prob,type="o",xlab="予備必要台数",main="更なる予備が必要になる確率",ylab="確率",col=2)

library(dplyr)

# kendall の結果 p.65

# 窓口が一つの場合で、到着分布を平均1/lambdaの指数分布、サービス分布を平均1/muの一般分布とする

N=100;ave=10

samples=rexp(N,rate=1/ave)

# 待ち時間の分布は指数分布なので最尤推定量で求める

lambda=1/mean(samples)

Data=data.frame(num=1:N,sample=samples,predict=rexp(N,rate=lambda))

plot(Data$num,Data$sample,type="p",col=2,xlim=c(min(Data[,1]),max(Data[,1])),ylim=c(min(Data[,2:3]),max(Data[,2:3])),xlab = "num",ylab="value")
par(new=T)
plot(Data$num,Data$predict,type="p",col=3,xlim=c(min(Data[,1]),max(Data[,1])),ylim=c(min(Data[,2:3]),max(Data[,2:3])),xlab = "num",ylab="value")


# サービス分布の平均と分散を指定する

ave=1;var=8;mu=1/ave

# n0は客が立ち去る前に並んでいた客の列の長さ、n1は以前から並んでいた客に対して、新規で到着したお客さまを含めたすべてのお客様の列の長さ

# p(n0=0)

if((lambda/mu) < 1){
  print("平衡状態")

prob=1-lambda/mu

print(paste0("p(n0=0)=",prob))

# 新規で到着するお客様の人数(r)の平均と分散を計算する

print(paste0("単位時間における窓口に到着する平均人数",lambda/mu))

print(paste0("単位時間における窓口に到着する人数の分散",(lambda/mu)+(lambda^2)*var))

# 系の長さの平均L=E[n0]

print(paste0("系の長さの平均L=",(lambda/mu)+((lambda/mu)^2+var*lambda^2)/(2*(1-lambda/mu))))

# W_q：待ち時間の平均

print(paste0("待ち時間の平均W_q＝",((lambda/mu)^2+var*lambda^2)/(2*lambda*(1-lambda/mu))))

}

library(dplyr)

# M/G/1 p.67

# B(v)が指数分布の場合。到着時間の分布を指数分布とする。

# 到着時間の平均時間aveからパラメータlambdaを計算

ave=1.1;lambda=1/ave;

# サービス時間の平均時間aveからパラメータmuを計算。分散varを計算する。

ave=1;mu=1/ave;var=(1/mu)^2;



prob_data=data.frame(wait_cus_num=0:100,prob=0)

for(j in 1:nrow(prob_data)){
  
prob_data$prob[j]=(mu/(lambda+mu))*(lambda/(lambda+mu))^prob_data$wait_cus_num[j]
 
}

plot(prob_data$wait_cus_num,prob_data$prob,type="o",col=2,xlab = "waiting customers",ylab="prob")


if((lambda/mu)<1){

    print(paste0("平衡状態：",lambda/mu))

  }

# prob_data$prob[prob_data$prob< 10^(-10)]=0

prob_data=prob_data %>% dplyr::filter(prob!=0)

# prob_data$prob[nrow(prob_data)]=prob_data$prob[nrow(prob_data)]+1-sum(prob_data$prob)

prob_length=length(prob_data$prob[prob_data$prob!=0])

P=matrix(0,nrow=prob_length, ncol=prob_length)

P[2,]=prob_data$prob;P[1,]=prob_data$prob;

for(j in 3:nrow(P)){
  
  vec=c()
  
  k=j-2
  
  vec=c(rep(0,k),prob_data$prob[1:(nrow(prob_data)-k)])
  
  P[j,]=vec
}





# 要素を個別に計算することで定常分布を保存する

stationary_prob=data.frame(generation=1:2000,vec=array(0,dim=c(1,prob_length)))

prob_vec=prob_data$prob

for(j in 1:nrow(stationary_prob)){
 
  for(l in 1:ncol(P)){
    
    stationary_prob[j,(l+1)]=sum(prob_vec*P[,l])
 
     }
  
  prob_vec=stationary_prob[j,2:ncol(stationary_prob)] 
  
}

stationary_num=stationary_prob$generation[abs(stationary_prob[,2]-(1-lambda/mu))==min(abs(stationary_prob[,2]-(1-lambda/mu)))]

stationary_probability=stationary_prob[stationary_num,2:ncol(stationary_prob)]

stationary_probability[length(stationary_probability)]=stationary_probability[length(stationary_probability)]+1-sum(stationary_probability)

# 定常分布が計算できたかどうか定常分布の条件を確認する（計算間違いの確認）

confirm_stationary=data.frame(num=0:100,stationary_prob=t(stationary_probability),confirm_prob=0)

colnames(confirm_stationary)=c("num","stationary_prob","confirm_prob")

confirm_stationary$confirm_prob[1]=1-lambda/mu

for(l in 1:(nrow(confirm_stationary)-2)){
  
  prob=stationary_probability[1]*P[1,(l+1)]
  
  for(k in 1:(l+1)){
    
  prob=prob+P[1,(l+2-k)]*stationary_probability[k+1]
    
  }
  
  confirm_stationary$confirm_prob[confirm_stationary$num==l]=prob
  
}

confirm_stationary$confirm_prob=floor(as.numeric(confirm_stationary$confirm_prob)*10^7)/(10^7)

confirm_stationary$stationary_prob=floor(as.numeric(confirm_stationary$stationary_prob)*10^7)/(10^7)

confirm_stationary=head(confirm_stationary,nrow(confirm_stationary)-2)

p=lambda/mu

confirm_stationary=confirm_stationary %>% dplyr::mutate(prob=(1-p)*p^(seq(0,(nrow(confirm_stationary)-1),1)))

# 第5章　待ち行列　G/G/1

# D/M/1 p.86~

# サービス率
mu=6;mu=1/mu

# 到着率
lam=9;lam=1/lam

k=1

# newton-method

f=function(x){

val=(1+x/(mu*k))^k-exp(x/lam)

return(val)

}

eta=10^(-5)

h=0.0001

z=-10

while(abs(f(z))>10^(-5))({

df=(f(z+h)-f(z))/h

z=z-eta*f(z)/df

print(f(z))

})

c=-(1-abs(z)/mu)

# 積分方程式の解

f=function(x){

val=1-(1-abs(z)/mu)*exp(-abs(z)*x)

return(val)

}

X=seq(0.01,1000,0.01)

plot(X,f(X))

# 5.2 El/Ek/1 p.87(確認中)

mu=3;mu=1/mu

lam=6;lam=1/lam

l=5;k=8

pthi=l*lam/(k*mu)

dim=rep(0,l+k+1)

dim[k+1]=-(1+pthi);dim[l+k+1]=1;dim[1]=pthi

z=polyroot(dim)

norm=c()

for(i in 1:length(z)){

norm=c(norm,sqrt(Re(z[i])^2+Im(z[i])^2))

}

u=z[round(norm*1000)/1000<1]

v=z[norm>1]

A=array(0,dim=c(k,k))

vec=rep(1,k)

for(i in 1:k){

A[i,]=vec

vec=vec*u

}

p10=(1-lam/mu)/prod(1-v)

c_vec=c()

for(i in 1:k){

A_sub=A

A_sub[,i]=(1+pthi)^c(0:(k-1))

c_vec=c(c_vec,prod(svd(A_sub)$d)*(p10*pthi/prod(svd(A)$d)))

}