差分

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

この比較画面へのリンク

両方とも前のリビジョン前のリビジョン
次のリビジョン		前のリビジョン
r:maintenance:condition_monitoring [2018/12/14 22:12] – [マルコフ決定過程とベルマン更新] watalur: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>
 +
 +
 +=== 劣化と保全 ===
  
 故障時間のデータしかなけば、寿命分布を推定し、時間取り替えやブロック取り替えを考えるしかない。 故障時間のデータしかなけば、寿命分布を推定し、時間取り替えやブロック取り替えを考えるしかない。
行 253: 行 271:
 1000000/(1+r)*(1+r) 1000000/(1+r)*(1+r)
 </code> </code>
 +
 +現在の資産が、将来に幾らになるかをグラフに表してみた。
 +
 +<code>
 +r = 0.10 ; plot(c(0:10),(1+r)^(0:10),type="l",col="black",ylab="Future Value of 1",xlab="Future (Year)")
 +r = 0.05 ; lines(c(0:10),(1+r)^(0:10),col="blue")
 +r = 0.01 ; lines(c(0:10),(1+r)^(0:10),col="brown")
 +legend("topleft",col=c("black","blue","brown"),lty=c(1,1,1),legend=c("10%","5%","1%"))
 +</code>
 +
 +当たり前のことだが、利回りが低ければ、遠い将来にもあまり価値の上昇は見られない。利回りが大きくなればなるほど、その上昇幅も大きくなる。
 +
 +{{:r:maintenance:future-values.png|}}
 +
 +
 +同じ費用が将来に発生する場合、現在価値に直すとどれぐらい減少するかも、グラフに表してみた。
 +
 +<code>
 +r = 0.10 ; plot(c(0:10),1/(1+r)^(0:10),type="l",col="black",ylab="Present Value of 1",xlab="Future (Year)")
 +r = 0.05 ; lines(c(0:10),1/(1+r)^(0:10),col="blue")
 +r = 0.01 ; lines(c(0:10),1/(1+r)^(0:10),col="brown")
 +legend("bottomleft",col=c("black","blue","brown"),lty=c(1,1,1),legend=c("10%","5%","1%"))
 +</code>
 +
 +利回りは10%と5%と1%の3種類の場合に、現在がt=0で10期先まで計算してある。利回りが10%の場合は、同じ費用が10期先に発生しても、現在価値に直すと半分以下になる。
 +
 +{{:r:maintenance:present-values.png|}}
  
 一期間が短く、また近い将来に発生する費用のみを考えれば良い場合には、利回りを考慮せずに、将来の費用をそのまま現在価値とすることもある。 一期間が短く、また近い将来に発生する費用のみを考えれば良い場合には、利回りを考慮せずに、将来の費用をそのまま現在価値とすることもある。
行 300: 行 345:
 </code> </code>
 と、毎期の将来5期間分だけの和を計算する。 と、毎期の将来5期間分だけの和を計算する。
-=== マルコフ決定過程とベルマン更新 ===+=== 現在価値の総費用の期待値に基づく意思決定 ===
  
-現在の状態と、将来の毎時点における状態と各行動の選択の帰結として、将来に渡って発生する費用を最小にするように、現在の行動を選択するモデルが、マルコフ決定過程である。+現在の状態と、将来の毎時点における状態と各行動の選択の帰結として、将来に渡って発生する費用を最小にするように、現在の行動を選択するモデルが、[[::markov_decision_process|マルコフ決定過程]]である。
  
 例えば、現在の状態が2とする。ここで稼働継続を選択すると、費用0が発生する。もし、次の時点で状態が3に推移するとして、稼働継続を選択していれば費用2000が生じてしまう。しかし、次の時点の状態は2のままかもしれない。あるいは現時点で、取替を選択すれば費用150が生じるが、次の時点では状態が1に戻るので稼働継続を選択できる。この時、稼働継続と修理と、どちらが望ましいだろうか。 例えば、現在の状態が2とする。ここで稼働継続を選択すると、費用0が発生する。もし、次の時点で状態が3に推移するとして、稼働継続を選択していれば費用2000が生じてしまう。しかし、次の時点の状態は2のままかもしれない。あるいは現時点で、取替を選択すれば費用150が生じるが、次の時点では状態が1に戻るので稼働継続を選択できる。この時、稼働継続と修理と、どちらが望ましいだろうか。
行 322: 行 367:
  
 これで、稼働継続の方が安いことが分かる。 これで、稼働継続の方が安いことが分かる。
-この計算を将来におけるすべての可能性をカバーするように行い、費用の現在価値での総和の期待値を最小にする行動を選択するのが、マルコフ決定過程に基づく最適な行動選択となる。+この計算を将来におけるすべての可能性をカバーするように行い、費用の現在価値での総和の期待値を最小にする行動を選択するのが、[[::markov_decision_process|マルコフ決定過程]]に基づく最適な行動選択となる。 
 + 
 +[[::markov_decision_process|マルコフ決定過程]]には次の3種類がある。 
 +  * 保全を繰り返して無限に運用を続ける際の将来に渡る総費用の期待値を最小化する問題 
 +  * 有限期間で運用が終わるため、通期の総費用を最小化する問題 
 +  * 固定の長さの直近の有限期間の総費用の期待値を最小化する問題 
 + 
 +また将来の費用を現在価値に割り戻すか否かで、次の２つの場合がある。 
 +  * 割引あり 
 +  * 割引なし 
 + 
 +割引なしで無限期間、と言う組み合わせのみ、総期待費用が♾に発散するため、合理的な意思決定ができない。 
 + 
 +==== 劣化データの等間隔の状態への変換 ==== 
 + 
 +以下のコードは、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>