Applied Statistics and Data Mining

Rの上で、次の9行を実行しておくこと。（学外のコンピュータでは、最初の3行は不要）

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)

故障時間のデータしかなけば、寿命分布を推定し、時間取り替えやブロック取り替えを考えるしかない。 しかし、多くの製品に多くの時点での劣化状態を測定したデータが得られていれば、それに基づいて劣化状態の変化のモデルを作り、時間取り替えやブロック取り替え以外の保全方式も検討できるようになる。

• National Bridge Inventoryは、アメリカ合衆国のすべての橋の棚卸しのデータである。毎年の春にすべての橋の状況が登録され、公開されている。橋の所在地、役割、種類、主な建材も記されているが、保全の研究者にとってはEXCELLENT, VERY GOOD, GOOD, SATISFACTORY, FAIR, POOR, SERIOUS, CRITICAL, “IMMINENT” FAILURE, FAILEDの10段階の状態のデータがとても興味深い。
• Hard Drive Data and Statsは、アメリカにあるBlackblazeというオンラインバックアップサービスやクラウドストレージサービスを展開する会社が公開している、ハードディスクの信頼性試験のデータである。実際い運用しているデータセンターでの統計情報のみでなく、個別のデータまで公開されているところが、画期的と言える。
• http://orbit.dtu.dk/files/7728815/ris_r_1235.pdf
• https://www.ntnu.edu/ross/info/data

対象の状態を監視し、状態に応じて適切な保全を施すことを、状態監視保全という。 英語では状態監視がcondition monitoring、監視した状態に基づく保全がcondition based maintenanceとなる。状態監視保全は「状態に基づく保全と状態監視」という意味のcondition based maintenance and monitoring、condition based maintenance and condition monitoring、あるいはcondition monitoring and condition based maintenanceなどが使われるが、単にcondition based maintenanceと呼ばれることも多い。

状態が正常と故障(機能喪失、性能低下を含む)の2状態のみの場合の保全行動には、次の3種類しかない。

1. 正常な状態の製品を正常な代替品で交換する予防取替
2. 故障した状態の製品を正常な代替品で交換する事後取替
3. 故障した状態の製品を正常な状態に戻す修理

取替は交換の言い換えで、修理は状態の回復を意味する。もし状態の種類が正常と故障の間をもっと細かく、例えば劣化しているが使える状態、という新たな状態も含めて3種類の観測ができるなら、状況はもっと複雑になる。上の3つの保全行動に加えて、次の4つが加わる。

1. 故障した状態の製品を劣化している状態に戻す修理
2. 劣化している製品を正常な状態に戻す予防修理
3. 故障した状態の製品を劣化している製品で交換する事後取替
4. 劣化している製品を正常な状態に戻す予防取替

更に複雑な製品や大規模な製品では、いきなり製品全体の取替と全体を使い続ける2択問題では自由度が低い。むしろ部分ごとに状態監視と修理や取替が可能となるように、コンポーネントと呼ばれる単位を定めて、システムとして製品全体を設計する。このとき状態も１次元ではなく、コンポーネントの数だけのベクトルになり、保全行動もコンポーネントごとに検討する。また一つを修理するついでに別のコンポーネントも修理するなど、保全行動の経済性も考慮するようになっていく。以下では、現実の問題よりは少し単純な状況を想定して、状態監視保全が有効な場合の保全行動の合理的な選択について考えてもらう。

マルコフ過程(以下では離散時間離散状態マルコフ過程, あるいはマルコフ連鎖を想定)は、未来の状態を予測するには現在の状態を知っていれば十分であり、過去の状態推移の履歴は必要ない、というマルコフ性を満たす。マルコフ性とは、時点tの状態X[t]の確率分布がt以前の過去の状態推移の履歴H[t-]={x[t-1], x[t-2], …}の条件付き確率分布F(X[t]|H[t-])で与えられているとして、その確率分布が直前の時点の状態x[t-1]のみに依存するF(X[t]|H[t-])=F(X[t]|X[t-1]=x[t-1])という性質である。

劣化データがあれば、保全の対象システムの劣化の進行がマルコフ性を満たしているかどうかを確認できる。マルコフ性を満たしていれば、マルコフ連鎖に基づく劣化の推移のモデル化が有効である。以下では定常マルコフ過程と呼ばれる、２時点間の状態推移を表す条件付き確率分布が時点によらず共通な場合のみを扱う。

修理も取替もせずに壊れるまで使用することと、システムの状態が１次元の数直線上の点で表せるとき、システムのマルコフ連鎖は次のような推移行列を持つ。

システムの状態が3つに分類され、1が正常、2が劣化、3が故障とする。 次の遷移行列は、毎時点、90%の確率で同じ状態に留まり10%の確率で状態が悪化すること、 そして故障すると何も働きかけない限りその状態(3)に留まること、を表す。

これを行列で表すと、次のようになる。

これを、劣化を見守ったときに i→j と状態が変化する確率を、配列の第(i,j)成分とする、という意味で、次のように記していく。

i→jという遷移が発生するための条件の一つは、状態がiにあること、もう一つは、保全行動を取らずに劣化の進行を許すこと、である。P(劣化)[i→j]は、この二つの条件の下での条件付き確率であり、長く書けば「P[i→j|現時点で状態がiであること & 保全行動を取らずに劣化の進行を許すこと]」となる。

これをRのコードで表すと

P.Dgr = matrix(c(
  9/10, 1/10, 0,
  0, 9/10, 1/10,
  0, 0, 1), nrow=3, ncol=3, byrow=TRUE)

参考までに、上のグラフは次のコードで描いた。

require(igraph)
g = graph.adjacency(P.Dgr*10,weighted=TRUE)
E(g)$weight = c(2,1,2,1,2)
plot(g,edge.width=E(g2)$weight,vertex.color="green",vertex.label.family="Helvetica")

壊れるまで使い続ける場合は、マルコフ解析が有効である。

上の推移行列は、劣化が徐々に進行することを表している。 これは、状態回復について何も働きかけない、という行動の下での遷移行列と言える。

取替という保全を行うと、どのような状態からでも正常な状態(1)に推移する。 これも遷移行列で表現できる。

システムの状態が3つに分類され、1が正常、2が劣化、3が故障とする。 次の遷移行列は、どの状態でも取替を実施すると必ず、次の時点では状態1になることを表す。

これを、必ず取替を行なうことと状態がiにあることの条件の下で状態が i→j と変化する確率を、配列の第(i,j)成分とする、という意味で、次のように記していく。

これをRのコードで表すと

P.Rpl = matrix(c(
  1, 0, 0,
  1, 0, 0,
  1, 0, 0), nrow=3, ncol=3, byrow=TRUE)

となる。ただし、この遷移行列はただでは使えない。取替費用を投じて始めて、推移行列を最初のP.Kから、このP.Rplに替えることができる。

同様に、修理による状態の変化もマルコフ連鎖の遷移行列で表現できる。

システムの状態が3つに分類され、1が正常、2が劣化、3が故障とする。 次の遷移行列は、どの状態でも取替を実施すると必ず、次の時点では状態が1だけ回復することを表す。

これを、必ず修理を行なうことと状態がiにあることの２つの条件の下で状態が i→j と変化する確率を、配列の第(i,j)成分とする、という意味で、次のように記していく。

これをRのコードで表すと

P.Rpr = matrix(c(
  1, 0, 0,
  1, 0, 0,
  0, 1, 0), nrow=3, ncol=3, byrow=TRUE)

となる。ただし、この遷移行列はただでは使えない。修理費用を投じて始めて、推移行列を最初のP.Dgrから、このP.Rprに替えることができる。取替と比べて、修理は状態回復の幅が小さく、通常は費用が取替よりも安く済む。

他に、1/2の確率で一つ状態をよくするが、1/2の確率でそのままの状態に留まる、確率修理もある。

P.Rpr.P = matrix(c(
  1, 0, 0,
  1/2, 1/2, 0,
  0, 1/2, 1/2), nrow=3, ncol=3, byrow=TRUE)

同様に、上のグラフは次のコードで描いている。

g <- graph.adjacency(P.Rpl,weighted=TRUE)
E(g)$weight = c(1,1,1)
plot(g,edge.width=E(g)$weight,vertex.color="green",vertex.label.family="Helvetica")
g <- graph.adjacency(P.Rpr,weighted=TRUE)
E(g)$weight = c(2,2,2)
plot(g,edge.width=E(g)$weight,vertex.color="green",vertex.label.family="Helvetica")
g <- graph.adjacency(P.Rpr.P*2,weighted=TRUE)
E(g)$weight = c(2,1,1,1,1)
plot(g,edge.width=E(g)$weight,vertex.color="green",vertex.label.family="Helvetica")

費用や損失は小さくすることが望ましく、報酬や利得は大きくすることが望ましい。 保全では通常、費用を考えることが多い。 費用の頭文字をとって「C」で表す。

典型的な保全行動の種類は次の通りである。

状態を１戻す修理(小修理)のに、2→1と3→2で同じ部品を交換するとも限らない。そのため、修理による状態推移2→1と別の修理3→2同じ費用がかかるとは限らない。

また状態を１に戻す修理(大修理)では、2→1と3→1で通常は前者の方が費用が安い。

取替の費用は、現在の状態のシステムを下取りに出したり、処分して上がる利益がなければ、一定となる。

例えば、次のようなコスト設定があるとする。

これをRでは、例えば次のようなコードで与える。

C.Opr = c(0,0,2000)
C.Rpr = c(10,50,250) 
C.Ovh = c(20,70,180)
C.Rpl = c(150,150,150)
Cost = cbind(C.Opr, C.Rpr, C.Ovh, C.Rpl)
colnames(Cost) = c("Keep","Repair","Overhaul","Replace")
rownames(Cost) = c("1","2","3")

  Keep  Repair Overhaul Replace
1     0     10       20     150
2     0     50       70     150
3   2000    250      180     150

グラフは次のコードで描ける。

matplot(Cost,type="b")
legend("topleft",lty=c(1,2,3,4),
       col=c(1,2,3,4),
       legend=c("1:Keep", "2:Repair", "3:Overhaul", "4: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万円の支出が見込まれている時、この支出に備えるにはいくらを用意すればいいか、またそれを運用した結果、いくらになるかは、次のように計算できる。

r = 0.01
1000000/(1+r)
1000000/(1+r)*(1+r)

現在の資産が、将来に幾らになるかをグラフに表してみた。

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%"))

当たり前のことだが、利回りが低ければ、遠い将来にもあまり価値の上昇は見られない。利回りが大きくなればなるほど、その上昇幅も大きくなる。

同じ費用が将来に発生する場合、現在価値に直すとどれぐらい減少するかも、グラフに表してみた。

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%"))

利回りは10%と5%と1%の3種類の場合に、現在がt=0で10期先まで計算してある。利回りが10%の場合は、同じ費用が10期先に発生しても、現在価値に直すと半分以下になる。

一期間が短く、また近い将来に発生する費用のみを考えれば良い場合には、利回りを考慮せずに、将来の費用をそのまま現在価値とすることもある。

例えば、最初の10期は稼働を見守り、11期目に修理をし、また次の10期に稼働を見守り、今度は22期目に取替をしたとする。これらの行動選択の結果、次のように状態が推移したとする。

これらの費用を、現在価値に直すと、それぞれ(1+r)^kで割れば良い。(k=1,…,22) それをRで行うと、次のようになる。

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)

これらの総和が、現在価値に換算した将来の22期間分の総期待費用である。

総和は関数sumで計算できる。

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))

更に、常に将来の5期間分だけの費用を考え続けるなら、

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)))
}

と、毎期の将来5期間分だけの和を計算する。将来が長く、選択肢も多く、意思決定のための計算量が膨大になる場合に、検討する将来の期間の長さをこのように制限することがある。

もし、現在価値に換算するのに利回りを用いなければ、単に発生する費用の合計となる。

sum(c(0,0,0,0,0,0,0,0,0,0,50,0,0,0,0,0,0,0,0,0,0,150))

利回りを考えない場合も、常に将来の5期間分だけの費用を考え続けるなら、

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)]))
}

と、毎期の将来5期間分だけの和を計算する。

現在の状態と、将来の毎時点における状態と各行動の選択の帰結として、将来に渡って発生する費用を最小にするように、現在の行動を選択するモデルが、マルコフ決定過程である。

例えば、現在の状態が2とする。ここで稼働継続を選択すると、費用0が発生する。もし、次の時点で状態が3に推移するとして、稼働継続を選択していれば費用2000が生じてしまう。しかし、次の時点の状態は2のままかもしれない。あるいは現時点で、取替を選択すれば費用150が生じるが、次の時点では状態が1に戻るので稼働継続を選択できる。この時、稼働継続と修理と、どちらが望ましいだろうか。

このことを考えるために、将来の総期待費用を考える。 状態2で稼働を継続した場合には、確率9/10で次の時点でも状態は2だが、確率1/10で次の時点の状態は3になる。取替を選択すれば、確率1で次の時点の状態は1となる。

稼働と取替のいずれを選択するかを、１期間先までに発生する費用を比較して、安く済む方を選択したい。 そのために次の式で総費用の現在価値を比較する。

0+(0*9/10+2000*1/10)/(1+r)
150+0/(1+r)

これで、稼働継続の方が安いことが分かる。 この計算を将来におけるすべての可能性をカバーするように行い、費用の現在価値での総和の期待値を最小にする行動を選択するのが、マルコフ決定過程に基づく最適な行動選択となる。

マルコフ決定過程には次の3種類がある。

• 保全を繰り返して無限に運用を続ける際の将来に渡る総費用の期待値を最小化する問題
• 有限期間で運用が終わるため、通期の総費用を最小化する問題
• 固定の長さの直近の有限期間の総費用の期待値を最小化する問題

また将来の費用を現在価値に割り戻すか否かで、次の２つの場合がある。

• 割引あり
• 割引なし

割引なしで無限期間、と言う組み合わせのみ、総期待費用が♾に発散するため、合理的な意思決定ができない。

以下のコードは、Dataというオブジェクトに対して、いろいろ操作をするように書いてる。 そのため貰ったデータをオブジェクトDataに代入する。

Data = X

まずデータdataから、時点数n.epochとアイテム数n.itemを取得する。

n.item = dim(Data)[1]
n.epoch = dim(Data)[2]

次に、グラフを描いて、データの様子を確認する。

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,])
}

別のグラフの描き方もある。

matplot(t(Data))

最小値と最大値を捜しておく。

min(Data)
max(Data)

マルコフ決定過程で保全方策を最適化するには、このグラフのデータを離散の状態の間の推移に変換する必要がある。等間隔の区間で離散化するために、最小値status.min、最大値status.max、区間幅status.widthを定める。次の数値は、例であり、各自で設定する必要がある。

status.min = 0.65
status.max = 1.05
status.width = 0.1

status.minとstatus.maxは、離散化するためのデータの範囲を定めるので、データの最小値と最大値を含むように、少し広く設定するといい。上のコードはstatus.minを小数点以下第一位で切り捨て、status.maxを小数点以下第一位で切り上げしている。

status.min = floor(min(Data)*10)/10
status.max = ceiling(max(Data)*10)/10
status.width = 0.2

これで区間数n.statusを求めておく。

n.status = ceiling((status.max-status.min)/status.width)

以上の数字より、各区間の境界を次のように定める。

status.breaks = status.min + c(0:n.status)*status.width

これを用いて、離散化を行う。

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)
  }
}

次のグラフを描くと、以上の手順から幾つかの状態を遷移しているデータに変換できたことが見て取れる。

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,])
}

データの最小値と最大値を確認しておく。

min(Data)
max(Data)

最小値と最大値を参考に、次のように区間を設定する。

status.breaks = c(0.4,0.5,0.7,0.9,1.1)

状態数は4となるが、これもRに数えさせておく。

n.status = length(status.breaks)-1

これを用いて、離散化をしてもよい。

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,])
}

マルコフ連鎖の遷移行列は、markovchainパッケージの中の関数markovchainFitで推定できる。

markovchainFit(Data.states)

この推定した遷移行列を、マルコフ決定過程で用いるので、オブジェクトに代入しておく。

Data.mc = markovchainFit(Data.states)

どのような数値や情報が出力されたかは、関数strで確認できる。

str(Data.mc)

推定された遷移行列は、次の１行で取り出せる。

Data.mc$estimate@transitionMatrix

これをオブジェクトに保管しておく。

P.Dgr = Data.mc$estimate@transitionMatrix
rownames(P.Dgr) = c(1:n.status)
colnames(P.Dgr) = c(1:n.status)

他の行動の遷移行列は次のように作る事ができる。

状態を1つだけ回復する修理は次のコードで生成できる。

n.translation = n.status-1
rbind(c(1,rep(0,n.status)),
            cbind(diag(rep(1,n.status)),0))

これを行と列のラベルを付けてオブジェクトP.Rprに保管する。

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)

どの状態でも、新品に戻す取替を表す遷移行列は次のコードで生成できる。

cbind(1,matrix(0,n.status,n.status-1))

これも行と列のラベルを付けてオブジェクトP.Rplに保管する。

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)

マルコフ決定過程の定義に必要なものは次の２つ。

• 行動ごとの遷移行列(推移確率行列)P
• 行動と状態ごとの費用行列R

Rは手作業での設定になる。

MDPtoolboxパッケージのために、行動ごとの遷移行列をすべて１つの配列に収める。 次のようにオブジェクトPを設定する。

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

P.Dgrを正しく推定できていれば、Pはn.status×n.status×n.actionの3次元配列になる。

同じく、MDPtoolboxパッケージのために費用関数を設定する。ここは、以下のコードをそのまま貼るのではなく、数字の数を状態数に合わせて適宜変更する必要がある。行動の数n.actionは3としている。また行動の順序は、上の遷移行列の配列と同じにする必要がある。

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

6状態であれば次のようにする。

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

各自の費用関数は、下記の数値の数を入れ替えて設定することになる。

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

C.Dgr, C.Rpr, C.Rplそれぞれの定義行のc()の中に数字がn.status個並び、数字と数字の間にはカンマがなければならない。

以下の２つのアルゴリズムの詳細はmarkov_decision_processに譲る。 ここでは割引係数(=1/(1+利回り))を0.9に設定して、最適方策を試算している。

マルコフ決定過程の最適方策を価値反復によって求めるには、次の一行を実行すればよい。

mdp_value_iteration(P, R, 0.9)

マルコフ決定過程の最適方策を方策反復によって求めるには、次の一行を実行すればよい。

mdp_policy_iteration(P, R, 0.9)