差分

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

この比較画面へのリンク

両方とも前のリビジョン前のリビジョン
次のリビジョン		前のリビジョン
r:maintenance:condition_monitoring [2018/12/14 21:13] – [保全の行動の費用もしくは報酬] 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>
 +
 +
 +=== 劣化と保全 ===
  
 故障時間のデータしかなけば、寿命分布を推定し、時間取り替えやブロック取り替えを考えるしかない。 故障時間のデータしかなけば、寿命分布を推定し、時間取り替えやブロック取り替えを考えるしかない。
行 208: 行 226:
 |1|0|10|20|150| |1|0|10|20|150|
 |2|0|50|70|150| |2|0|50|70|150|
-|3|500|250|180|150|+|3|2000|250|180|150|
  
 これをRでは、例えば次のようなコードで与える。 これをRでは、例えば次のようなコードで与える。
  
 <code> <code>
-C.Opr = c(0,0,500)+C.Opr = c(0,0,2000)
 C.Rpr = c(10,50,250)  C.Rpr = c(10,50,250) 
 C.Ovh = c(20,70,180) C.Ovh = c(20,70,180)
 C.Rpl = c(150,150,150) C.Rpl = c(150,150,150)
 Cost = cbind(C.Opr, C.Rpr, C.Ovh, C.Rpl) Cost = cbind(C.Opr, C.Rpr, C.Ovh, C.Rpl)
-colnames(Cost) = c("Oper.","Repair","Overhaul","Replace")+colnames(Cost) = c("Keep","Repair","Overhaul","Replace")
 rownames(Cost) = c("1","2","3") rownames(Cost) = c("1","2","3")
 </code> </code>
  
 <code> <code>
-  Oper. Repair Overhaul Replace+  Keep  Repair Overhaul Replace
 1     0     10       20     150 1     0     10       20     150
 2     0     50       70     150 2     0     50       70     150
-3   500    250      180     150+3   2000    250      180     150
 </code> </code>
  
行 242: 行 260:
 グラフからは、状態が1の場合には、状態に介入しない行動Keepの費用が一番安いことが見て取れる。このことから、状態1ではKeepを選択するのが望ましいように思われる。状態が2の場合には状態に介入しない行動Keepの費用が一番安いが、将来のことを考えると状態1に戻す行動Repairも気になる。状態が3の場合には状態1に戻す行動Replaceの費用が一番安いことが見て取れる。もし、状態3から状態2に戻す行動Repairの費用が、Replaceの費用よりも安ければ、どちらにしようか悩む。 グラフからは、状態が1の場合には、状態に介入しない行動Keepの費用が一番安いことが見て取れる。このことから、状態1ではKeepを選択するのが望ましいように思われる。状態が2の場合には状態に介入しない行動Keepの費用が一番安いが、将来のことを考えると状態1に戻す行動Repairも気になる。状態が3の場合には状態1に戻す行動Replaceの費用が一番安いことが見て取れる。もし、状態3から状態2に戻す行動Repairの費用が、Replaceの費用よりも安ければ、どちらにしようか悩む。
  
-この意思決定問題の解を、将来に生じる費用の合計に基づいて合理的に選択するように得るためのモデルが、マルコフ決定過程である。+この意思決定問題の解を、将来に生じる費用の合計に基づいて合理的に選択するように得るためのモデルが、マルコフ決定過程である。マルコフ決定過程は複数ある行動の選択が、対象にマルコフ行列で表される影響を与える場合に、総期待費用に基づいて状態に応じた意思決定を行う場合に有効なモデルである。
  
-=== マルコフ決定過程 ===+=== 現在価値 ===
  
-複数ある行動を、状態に応じて施すことを考える場合は、マルコフ決定過程が有効である。+将来に生じる費用を現在の支出と比較するために、その支出を現在から備えておくとしたらどれぐらいの金額が必要か、という金額に換算した金額を現在価値という。1期間あたり利回りがrの資産運用が可能な状況では、1期間後に発生する費用Cの現在価値は C/(1+r) となる。例えば利回りをr=0.01、一期先に100万円の支出が見込まれている時、この支出に備えるにはいくらを用意すればいいか、またそれを運用した結果、いくらになるかは、次のように計算できる。 
 + 
 +<code> 
 +r = 0.01 
 +1000000/(1+r) 
 +1000000/(1+r)*(1+r) 
 +</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|}} 
 + 
 +一期間が短く、また近い将来に発生する費用のみを考えれば良い場合には、利回りを考慮せずに、将来の費用をそのまま現在価値とすることもある。 
 + 
 +=== 総期待費用 === 
 + 
 +例えば、最初の10期は稼働を見守り、11期目に修理をし、また次の10期に稼働を見守り、今度は22期目に取替をしたとする。これらの行動選択の結果、次のように状態が推移したとする。 
 + 
 +|k|1|2|3|4|5|6|7|8|9|10|11|12|13|14|15|16|17|18|19|20|21|22| 
 +|行動|K|K|K|K|K|K|K|K|K|K|Rpr|K|K|K|K|K|K|K|K|K|K|Rpl| 
 +|状態|1|1|1|1|1|1|1|2|2|2|1|1|1|1|1|1|1|1|1|2|2|1| 
 +|費用|0|0|0|0|0|0|0|0|0|0|50|0|0|0|0|0|0|0|0|0|0|150| 
 + 
 +これらの費用を、現在価値に直すと、それぞれ(1+r)^kで割れば良い。(k=1,...,22) 
 +それをRで行うと、次のようになる。 
 +<code> 
 +r = 0.01 
 +c(0,0,0,0,0,0,0,0,0,0,50,0,0,0,0,0,0,0,0,0,0,150)/(1+r)^c(1:22) 
 +</code> 
 +これらの総和が、現在価値に換算した将来の22期間分の総期待費用である。 
 + 
 +総和は関数sumで計算できる。 
 + 
 +<code> 
 +sum(c(0,0,0,0,0,0,0,0,0,0,50,0,0,0,0,0,0,0,0,0,0,150)/(1+r)^c(1:22)) 
 +</code> 
 + 
 +更に、常に将来の5期間分だけの費用を考え続けるなら、 
 +<code> 
 +for( i in c(1:18) ) { 
 +  print(sum(c(0,0,0,0,0,0,0,0,0,0,50,0,0,0,0,0,0,0,0,0,0,150)[i:(i+4)]/(1+r)^c(1:5))) 
 +} 
 +</code> 
 +と、毎期の将来5期間分だけの和を計算する。将来が長く、選択肢も多く、意思決定のための計算量が膨大になる場合に、検討する将来の期間の長さをこのように制限することがある。 
 + 
 +もし、現在価値に換算するのに利回りを用いなければ、単に発生する費用の合計となる。 
 + 
 +<code> 
 +sum(c(0,0,0,0,0,0,0,0,0,0,50,0,0,0,0,0,0,0,0,0,0,150)) 
 +</code> 
 + 
 +利回りを考えない場合も、常に将来の5期間分だけの費用を考え続けるなら、 
 +<code> 
 +for( i in c(1:18) ) { 
 +  print(sum(c(0,0,0,0,0,0,0,0,0,0,50,0,0,0,0,0,0,0,0,0,0,150)[i:(i+4)])) 
 +} 
 +</code> 
 +と、毎期の将来5期間分だけの和を計算する。 
 +=== 現在価値の総費用の期待値に基づく意思決定 === 
 + 
 +現在の状態と、将来の毎時点における状態と各行動の選択の帰結として、将来に渡って発生する費用を最小にするように、現在の行動を選択するモデルが、[[::markov_decision_process|マルコフ決定過程]]である。 
 + 
 +例えば、現在の状態が2とする。ここで稼働継続を選択すると、費用0が発生する。もし、次の時点で状態が3に推移するとして、稼働継続を選択していれば費用2000が生じてしまう。しかし、次の時点の状態は2のままかもしれない。あるいは現時点で、取替を選択すれば費用150が生じるが、次の時点では状態が1に戻るので稼働継続を選択できる。この時、稼働継続と修理と、どちらが望ましいだろうか。 
 + 
 +このことを考えるために、将来の総期待費用を考える。 
 +状態2で稼働を継続した場合には、確率9/10で次の時点でも状態は2だが、確率1/10で次の時点の状態は3になる。取替を選択すれば、確率1で次の時点の状態は1となる。 
 + 
 +|行動選択|現時点にかかる費用|１時点先の状態の候補|１時点先の行動選択の候補|１期先までの費用の期待値の合計| 
 +|稼働|0|2 または 3|2なら稼働、3なら取替|0+(0*9/10+2000*1/10)/(1+r)| 
 +|取替|150|1|稼働|150+0/(1+r)| 
 + 
 +稼働と取替のいずれを選択するかを、１期間先までに発生する費用を比較して、安く済む方を選択したい。 
 +そのために次の式で総費用の現在価値を比較する。 
 + 
 +<code> 
 +0+(0*9/10+2000*1/10)/(1+r) 
 +150+0/(1+r) 
 +</code> 
 + 
 +これで、稼働継続の方が安いことが分かる。 
 +この計算を将来におけるすべての可能性をカバーするように行い、費用の現在価値での総和の期待値を最小にする行動を選択するのが、[[::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>