文書の過去の版を表示しています。
+目次
時間計画保全
寿命分布の設定
午前の結果に基づいて分布を選択する。 対応する確率分布のpdensfun, ddensfun, rdensfunをRに貼り付ける。
# 寿命分布の累積分布関数 pdensfun = pweibull # ワイブル分布 pdensfun = plnorm # 対数正規分布 pdensfun = pgamma # ガンマ分布 # 寿命分布の確率密度関数 ddensfun = dweibull # ワイブル分布 ddensfun = dlnorm # 対数正規分布 ddensfun = dgamma # ガンマ分布 # 寿命分布の乱数 rdensfun = rweibull # ワイブル分布 rdensfun = rlnorm # 対数正規分布 rdensfun = rgamma # ガンマ分布
パラメータの推定値を与える。
# 寿命分布のパラメータ 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)
確率分布の密度関数と生存関数
コストレートの式の分子の被積分関数が確率密度関数なので、パラメータまで与えた関数fを定義しておく。
# 確率密度関数 f <- function(x) { return(ddensfun(x,rdensfun.params[1],rdensfun.params[2])) }
コストレートの式の分母の被積分関数が生存関数(=1-累積分布関数のこと)なので、これもパラメータまで与えた関数F.barを定義しておく。
# 生存関数 F.bar <- function(x) { return(1-pdensfun(x,rdensfun.params[1],rdensfun.params[2])) }
時間取り替え
コストレートは、分母も分子も積分を含んでいる。上に紹介した関数integrateと、被積分関数の与え方に従えば、次のようにコストレートをRの関数で表現できる。
# コストレートの式 (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)) }
例えば
g(1)
と入力すると、T=1の場合のコストレートが計算される。
コストレートの最小化に、関数nlminbを使うことができる。nlminbは最小化をしてくれるRの関数で、初期値T=1から始めて、最小値を探すには、次の1行を実行する。
nlminb(1,g)
最小解のみでなく、どのような関数を最小化したかを確認するには、グラフも描いてみると良いだろう。
plot(g,xlim=c(0,20))
他に、次のようにしても、同じグラフを描くことができる。
# コストレートのグラフを描く g.list <- NULL for( i in c(1:200) ) { g.list <- append(g.list,g(i/10)) } plot(c(1:500)/10,g.list,type="b")
0.1ずつ関数の値を評価して、折れ線グラフを描いている。
ブロック取り替え
# 再生関数をモンテカルロ近似で求める # 再生関数が数式で表せる分布の方が少ない # 準備のための乱数データ発生 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))