差分

このページの2つのバージョン間の差分を表示します。

この比較画面にリンクする

両方とも前のリビジョン 前のリビジョン
次のリビジョン 		前のリビジョン
r:maintenance:condition_monitoring [2018/12/16 16:39]
watalu [現在価値の総費用の期待値に基づく意思決定] 		r:maintenance:condition_monitoring [2019/12/16 16:57] (現在)
watalu [費用関数の設定]
行 1: 行 1:
 ==== 状態監視保全 ==== ==== 状態監視保全 ====
 +
 +=== 準備 ===
 +
 +Rの上で、次の9行を実行しておくこと。（学外のコンピュータでは、最初の3行は不要）
 +<code>
 +Sys.setenv("http_proxy"="http://130.153.8.19:8080/")
 +Sys.setenv("ftp_proxy"="http://130.153.8.19:8080/")
 +Sys.setenv("https_proxy"="http://130.153.8.19:8080/")
 +install.packages("markovchain")
 +install.packages("igraph")
 +install.packages("MDPtoolbox")
 +library(markovchain)
 +library(igraph)
 +library(MDPtoolbox)
 +</code>
 +
 +
 +=== 劣化と保全 ===
  
 故障時間のデータしかなけば、寿命分布を推定し、時間取り替えやブロック取り替えを考えるしかない。 故障時間のデータしかなけば、寿命分布を推定し、時間取り替えやブロック取り替えを考えるしかない。
行 361: 行 379:
  
 割引なしで無限期間、と言う組み合わせのみ、総期待費用が♾に発散するため、合理的な意思決定ができない。 割引なしで無限期間、と言う組み合わせのみ、総期待費用が♾に発散するため、合理的な意思決定ができない。
 +
 +==== 劣化データの等間隔の状態への変換 ====
 +
 +以下のコードは、Dataというオブジェクトに対して、いろいろ操作をするように書いてる。
 +そのため貰ったデータをオブジェクトDataに代入する。
 +
 +<code>
 +Data = X
 +</code>
 +
 +まずデータdataから、時点数n.epochとアイテム数n.itemを取得する。
 +
 +<code>
 +n.item = dim(Data)[1]
 +n.epoch = dim(Data)[2]
 +</code>
 +
 +次に、グラフを描いて、データの様子を確認する。
 +
 +<code>
 +plot(c(1,n.epoch),c(min(Data),max(Data)),
 +     xlab="Time",ylab="Data",type="n")
 +for( i in c(1:n.item) ) {
 +  lines(c(1:n.epoch),Data[i,])
 +}
 +</code>
 +
 +別のグラフの描き方もある。
 +
 +<code>
 +matplot(t(Data))
 +</code>
 +
 +最小値と最大値を捜しておく。
 +
 +<code>
 +min(Data)
 +max(Data)
 +</code>
 +
 +マルコフ決定過程で保全方策を最適化するには、このグラフのデータを離散の状態の間の推移に変換する必要がある。等間隔の区間で離散化するために、最小値status.min、最大値status.max、区間幅status.widthを定める。次の数値は、例であり、各自で設定する必要がある。
 +
 +<code>
 +status.min = 0.65
 +status.max = 1.05
 +status.width = 0.1
 +</code>
 +
 +status.minとstatus.maxは、離散化するためのデータの範囲を定めるので、データの最小値と最大値を含むように、少し広く設定するといい。上のコードはstatus.minを小数点以下第一位で切り捨て、status.maxを小数点以下第一位で切り上げしている。
 +
 +<code>
 +status.min = floor(min(Data)*10)/10
 +status.max = ceiling(max(Data)*10)/10
 +status.width = 0.2
 +</code>
 +
 +これで区間数n.statusを求めておく。
 +
 +<code>
 +n.status = ceiling((status.max-status.min)/status.width)
 +</code>
 +
 +以上の数字より、各区間の境界を次のように定める。
 +
 +<code>
 +status.breaks = status.min + c(0:n.status)*status.width
 +</code>
 +
 +これを用いて、離散化を行う。
 +
 +<code>
 +Data.states = Data
 +for( i in c(1:n.status) ) {
 +  if(sum((Data>status.breaks[i]) & (Data<=status.breaks[i+1]))>0) {
 +    Data.states[ (Data>status.breaks[i]) & (Data<=status.breaks[i+1]) ] = as.character(n.status-i+1)
 +  }
 +}
 +</code>
 +
 +次のグラフを描くと、以上の手順から幾つかの状態を遷移しているデータに変換できたことが見て取れる。
 +
 +<code>
 +plot(c(1,n.epoch),c(min(Data.states),max(Data.states)),
 +     xlab="Time",ylab="Status",type="n")
 +for( i in c(1:n.item) ) {
 +  lines(c(1:n.epoch),Data.states[i,])
 +}
 +</code>
 +
 +==== 劣化データの任意の間隔の状態への変換 ====
 +
 +データの最小値と最大値を確認しておく。
 +
 +<code>
 +min(Data)
 +max(Data)
 +</code>
 +
 +最小値と最大値を参考に、次のように区間を設定する。
 +
 +<code>
 +status.breaks = c(0.4,0.5,0.7,0.9,1.1)
 +</code>
 +
 +状態数は4となるが、これもRに数えさせておく。
 +
 +<code>
 +n.status = length(status.breaks)-1
 +</code>
 +
 +これを用いて、離散化をしてもよい。
 +
 +<code>
 +Data.states = Data
 +for( i in c(1:n.status) ) {
 +  if(sum((Data>status.breaks[i]) & (Data<=status.breaks[i+1]))>0) {
 +    Data.states[ (Data>status.breaks[i]) & (Data<=status.breaks[i+1]) ] = as.character(ceiling(n.status-i+1))
 +  }
 +}
 +plot(c(1,n.epoch),c(min(Data.states),max(Data.states)),
 +     xlab="Time",ylab="Status",type="n")
 +for( i in c(1:n.item) ) {
 +  lines(c(1:n.epoch),Data.states[i,])
 +}
 +</code>
 +
 +==== マルコフ連鎖の遷移行列の推定 ====
 +
 +マルコフ連鎖の遷移行列は、markovchainパッケージの中の関数markovchainFitで推定できる。
 +<code>
 +markovchainFit(Data.states)
 +</code>
 +
 +この推定した遷移行列を、マルコフ決定過程で用いるので、オブジェクトに代入しておく。
 +
 +<code>
 +Data.mc = markovchainFit(Data.states)
 +</code>
 +
 +どのような数値や情報が出力されたかは、関数strで確認できる。
 +
 +<code>
 +str(Data.mc)
 +</code>
 +
 +推定された遷移行列は、次の１行で取り出せる。
 +<code>
 +Data.mc$estimate@transitionMatrix
 +</code>
 +
 +これをオブジェクトに保管しておく。
 +
 +<code>
 +P.Dgr = Data.mc$estimate@transitionMatrix
 +rownames(P.Dgr) = c(1:n.status)
 +colnames(P.Dgr) = c(1:n.status)
 +</code>
 +
 +他の行動の遷移行列は次のように作る事ができる。
 +
 +状態を1つだけ回復する修理は次のコードで生成できる。
 +
 +<code>
 +n.translation = n.status-1
 +rbind(c(1,rep(0,n.status)),
 +            cbind(diag(rep(1,n.status)),0))
 +</code>
 +
 +これを行と列のラベルを付けてオブジェクトP.Rprに保管する。
 +<code>
 +P.Rpr = rbind(c(1,rep(0,n.status-1)),
 +            cbind(diag(rep(1,n.status-1)),0))
 +rownames(P.Rpr) = c(1:n.status)
 +colnames(P.Rpr) = c(1:n.status)
 +</code>
 +
 +どの状態でも、新品に戻す取替を表す遷移行列は次のコードで生成できる。
 +
 +<code>
 +cbind(1,matrix(0,n.status,n.status-1))
 +</code>
 +
 +これも行と列のラベルを付けてオブジェクトP.Rplに保管する。
 +
 +<code>
 +P.Rpl = cbind(1,matrix(0,n.status,n.status-1))
 +rownames(P.Rpl) = c(1:n.status)
 +colnames(P.Rpl) = c(1:n.status)
 +</code>
 +
 +
 +
 +==== MDPtoobox利用の準備 ====
 +
 +マルコフ決定過程の定義に必要なものは次の２つ。
 +
 +  * 行動ごとの遷移行列(推移確率行列)P
 +  * 行動と状態ごとの費用行列R
 +
 +Rは手作業での設定になる。
 +
 +===== 遷移行列の配列の設定 =====
 +
 +MDPtoolboxパッケージのために、行動ごとの遷移行列をすべて１つの配列に収める。
 +次のようにオブジェクトPを設定する。
 +
 +<code>
 +n.action = 3
 +P = array(dim=c(n.status, n.status, n.action))
 +P[,,1] = P.Dgr
 +P[,,2] = P.Rpr
 +P[,,3] = P.Rpl
 +</code>
 +
 +P.Dgrを正しく推定できていれば、Pはn.status×n.status×n.actionの3次元配列になる。
 +
 +===== 費用関数の設定 =====
 +
 +同じく、MDPtoolboxパッケージのために費用関数を設定する。ここは、以下のコードをそのまま貼るのではなく、数字の数を状態数に合わせて適宜変更する必要がある。行動の数n.actionは3としている。また行動の順序は、上の遷移行列の配列と同じにする必要がある。
 +
 +<code>
 +C.Dgr = c(0,0,2000)
 +C.Rpr = c(10,50,250) 
 +C.Rpl = c(150,150,150)
 +Cost = cbind(C.Dgr, C.Rpr, C.Rpl)
 +colnames(Cost) = c("Keep","Repair","Replace")
 +rownames(Cost) = c("1","2","3")
 +R = -Cost
 +</code>
 +
 +6状態であれば次のようにする。
 +
 +<code>
 +C.Dgr = c(0,0,0,0,0,2000)
 +C.Rpr = c(10,50,100,160,230,310) 
 +C.Rpl = c(150,150,150,150,150,150)
 +Cost = cbind(C.Dgr, C.Rpr, C.Rpl)
 +colnames(Cost) = c("Keep","Repair","Replace")
 +rownames(Cost) = as.character(1:6)
 +R = -Cost
 +</code>
 +
 +各自の費用関数は、下記の数値の数を入れ替えて設定することになる。
 +
 +<code>
 +C.Dgr = c()
 +C.Rpr = c() 
 +C.Rpl = c()
 +Cost = cbind(C.Dgr, C.Rpr, C.Rpl)
 +colnames(Cost) = c("Keep","Repair","Replace")
 +rownames(Cost) = as.character(1:n.status)
 +R = -Cost
 +</code>
 +
 +C.Dgr, C.Rpr, C.Rplそれぞれの定義行のc()の中に数字がn.status個並び、数字と数字の間にはカンマがなければならない。
 +==== 最適方策の算出 ====
 +
 +以下の２つのアルゴリズムの詳細は[[::markov_decision_process]]に譲る。
 +ここでは割引係数(=1/(1+利回り))を0.9に設定して、最適方策を試算している。
 +=== 価値反復法 ===
 +
 +マルコフ決定過程の最適方策を価値反復によって求めるには、次の一行を実行すればよい。
 +
 +<code>
 +mdp_value_iteration(P, R, 0.9)
 +</code>
 +
 +=== 方策反復法 ===
 +
 +マルコフ決定過程の最適方策を方策反復によって求めるには、次の一行を実行すればよい。
 +
 +<code>
 +mdp_policy_iteration(P, R, 0.9)
 +</code>
 +
 +