文書の過去の版を表示しています。
時間計画保全
寿命分布の設定
午前の結果に基づいて分布を選択し、パラメータの推定値を与える。
# 寿命分布 pdensfun = pweibull # ワイブル分布 pdensfun = plnorm # 対数正規分布 pdensfun = pgamma # ガンマ分布 # 寿命分布のパラメータ rdensfun.params = c(2.2357389,9.7422387)
コストの設定
#今回の想定パラメータ Cc = 7800+40000 Cp = 7800+5000
数値積分
時間取替のコストレートの式には、積分が含まれている。Rでは、関数integrateによる積分の機能が提供されている。
三角関数sinを0から2πまで積分してみる。
integrate(sin,0,2*pi)
Rの中で定義済みの関数は、このようにそのまま定積分を計算できる。中で定義されていない関数、あるいは積分変数以外の変数を持つ関数は、積分変数のみを引数として、被積分関数の値を返すRの関数を定義して、integrateに渡す。
f = function(x) { f.ret = dlnorm(x, meanlog=2, sdlog=1) }
このように定義した関数であれば、次のように範囲を指定して、グラフを描くことができる。
plot(f, xlim=c(0,20))
この関数を(0,♾)の範囲で定積分を行うには、次のようにRに指示すれば良い。
integrate(f,0,Inf)
生存関数
コストレートの式の分子の被積分関数が生存関数(=1-累積分布関数のこと)なので、その関数を定義しておく。
# 生存関数 F.bar <- function(x) { return(1-pdensfun(x,rdensfun.params[1],rdensfun.params[2])) }
時間取り替え
コストレートは、積分
# コストレートの式 (integrateは1変数関数の数値積分をしてくれるRの関数) g <- function(x) { return((Cc*integrate(f,0,x)$value+Cp*(1-integrate(f,0,x)$value))/(integrate(F.bar,0,x)$value)) }
# nlminbは最小化をしてくれるRの関数 nlminb(1,g)
# コストレートのグラフを描く g.list <- NULL for( i in c(1:500) ) { g.list <- append(g.list,g(i/10)) } plot(c(1:500)/10,g.list,type="b")
ブロック取り替え
# 再生関数をモンテカルロ近似で求める # 再生関数が数式で表せる分布の方が少ない # 準備のための乱数データ発生 M.t.prep <- function() { rdist <- function(n) { return(rdensfun(n,rdensfun.params[1],rdensfun.params[2])) } n <- 100000 X <- rep(0,n) K <- 50 F.k <- array(0,dim=c(K,n)) for( k in c(1:K) ) { X <- X + rdist(n) F.k[k,] <- X } return(F.k) } # 再生関数をモンテカルロ計算で求める M.t <- function(block.cycles) { n <- 100000 X <- rep(0,n) K <- 50 M.ret <- rep(0,length(block.cycles)) for( i in c(1:length(block.cycles)) ) { for( k in c(1:K) ) { M.ret[i] <- M.ret[i] + sum(M.t.data[k,]<=block.cycles[i])/n } } return(M.ret) }
# 使う前に一度だけ準備を行う(分布を変更する都度、実行しなおす必要あり) M.t.data = M.t.prep() # 再生関数のグラフ plot(1:100,M.t(1:100)) # ブロック取り替えのコストレートのグラフ plot(1:100,(M.t(1:100)*Cc+Cp)/(1:100), xlab="T", ylab="Cost Rate") # ブロック取り替えの子ストレート Cost.rate.BR <- function(cycle.length) { return((M.t(cycle.length)*Cc+Cp)/cycle.length) } # ブロック取り替えの最適化をRの関数optimizeで実行する optimize(Cost.rate.BR,c(0,10))