差分

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

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

両方とも前のリビジョン 前のリビジョン
次のリビジョン 		前のリビジョン
r:maintenance:condition_monitoring [2018/12/14 12:02]
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>
 +
 +
 +=== 劣化と保全 ===
  
 故障時間のデータしかなけば、寿命分布を推定し、時間取り替えやブロック取り替えを考えるしかない。 故障時間のデータしかなけば、寿命分布を推定し、時間取り替えやブロック取り替えを考えるしかない。
行 27: 行 45:
 更に複雑な製品や大規模な製品では、いきなり製品全体の取替と全体を使い続ける2択問題では自由度が低い。むしろ部分ごとに状態監視と修理や取替が可能となるように、コンポーネントと呼ばれる単位を定めて、システムとして製品全体を設計する。このとき状態も１次元ではなく、コンポーネントの数だけのベクトルになり、保全行動もコンポーネントごとに検討する。また一つを修理するついでに別のコンポーネントも修理するなど、保全行動の経済性も考慮するようになっていく。以下では、現実の問題よりは少し単純な状況を想定して、状態監視保全が有効な場合の保全行動の合理的な選択について考えてもらう。 更に複雑な製品や大規模な製品では、いきなり製品全体の取替と全体を使い続ける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])という性質である。 マルコフ過程(以下では離散時間離散状態マルコフ過程, あるいはマルコフ連鎖を想定)は、未来の状態を予測するには現在の状態を知っていれば十分であり、過去の状態推移の履歴は必要ない、というマルコフ性を満たす。マルコフ性とは、時点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])という性質である。
行 43: 行 61:
 次の遷移行列は、毎時点、90%の確率で同じ状態に留まり10%の確率で状態が悪化すること、 次の遷移行列は、毎時点、90%の確率で同じ状態に留まり10%の確率で状態が悪化すること、
 そして故障すると何も働きかけない限りその状態(3)に留まること、を表す。 そして故障すると何も働きかけない限りその状態(3)に留まること、を表す。
 +
 +{{:r:maintenance:trantisionmatrix-keep.png|}}
 +
 +これを行列で表すと、次のようになる。
  
 |P[1->1]=9/10|P[1->2]=1/10|P[1->3]=0| |P[1->1]=9/10|P[1->2]=1/10|P[1->3]=0|
行 63: 行 85:
   0, 0, 1), nrow=3, ncol=3, byrow=TRUE)   0, 0, 1), nrow=3, ncol=3, byrow=TRUE)
 </code> </code>
 +
 +参考までに、上のグラフは次のコードで描いた。
 +
 +<code>
 +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")
 +</code>
 +
 +=== マルコフ連鎖 ===
 +
 +壊れるまで使い続ける場合は、マルコフ解析が有効である。
 +
 +==== 保全行動の選択問題 ====
 === 保全行動の遷移行列による表現 === === 保全行動の遷移行列による表現 ===
  
行 82: 行 119:
 |P[2->1]=1|P[2->2]=0|P[2->3]=0| |P[2->1]=1|P[2->2]=0|P[2->3]=0|
 |P[3->1]=1|P[3->2]=0|P[3->3]=0| |P[3->1]=1|P[3->2]=0|P[3->3]=0|
 +
 +{{:r:maintenance:transitionmatrix-replace.png|}}
  
 これを、必ず取替を行なうことと状態がiにあることの条件の下で状態が i->j と変化する確率を、配列の第(i,j)成分とする、という意味で、次のように記していく。 これを、必ず取替を行なうことと状態がiにあることの条件の下で状態が i->j と変化する確率を、配列の第(i,j)成分とする、という意味で、次のように記していく。
行 111: 行 150:
 |P[2->1]=1|P[2->2]=0|P[2->3]=0| |P[2->1]=1|P[2->2]=0|P[2->3]=0|
 |P[3->1]=0|P[3->2]=1|P[3->3]=0| |P[3->1]=0|P[3->2]=1|P[3->3]=0|
 +
 +{{:r:maintenance:transitionmatrix-repair.png|}}
  
 これを、必ず修理を行なうことと状態がiにあることの２つの条件の下で状態が i->j と変化する確率を、配列の第(i,j)成分とする、という意味で、次のように記していく。 これを、必ず修理を行なうことと状態がiにあることの２つの条件の下で状態が i->j と変化する確率を、配列の第(i,j)成分とする、という意味で、次のように記していく。
行 126: 行 167:
 </code> </code>
 となる。ただし、この遷移行列はただでは使えない。修理費用を投じて始めて、推移行列を最初のP.Dgrから、このP.Rprに替えることができる。取替と比べて、修理は状態回復の幅が小さく、通常は費用が取替よりも安く済む。 となる。ただし、この遷移行列はただでは使えない。修理費用を投じて始めて、推移行列を最初のP.Dgrから、このP.Rprに替えることができる。取替と比べて、修理は状態回復の幅が小さく、通常は費用が取替よりも安く済む。
 +
 +他に、1/2の確率で一つ状態をよくするが、1/2の確率でそのままの状態に留まる、確率修理もある。
 +
 +{{:r:maintenance:transitionmatrix-repair-with-probability.png|}}
 +
 +<code>
 +P.Rpr.P = matrix(c(
 +  1, 0, 0,
 +  1/2, 1/2, 0,
 +  0, 1/2, 1/2), nrow=3, ncol=3, byrow=TRUE)
 +</code>
 +
 +同様に、上のグラフは次のコードで描いている。
 +
 +<code>
 +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")
 +</code>
 +
 +=== 保全の行動の費用もしくは報酬 ===
 +
 +費用や損失は小さくすることが望ましく、報酬や利得は大きくすることが望ましい。
 +保全では通常、費用を考えることが多い。
 +費用の頭文字をとって「C」で表す。
 +
 +典型的な保全行動の種類は次の通りである。
 +
 +|保全行動|説明|システムの状態への影響行使|
 +|稼働(継続)|稼働を見守る|状態遷移に影響を行使せず、劣化過程に委ねる|
 +|確率修理|点検・診断の後、状態に応じて定められた修理を行う|劣化過程を止めて、修理内容に応じた状態に戻すが、戻る先は幾つかの状態のいずれかとなる|
 +|小修理|点検・診断の後、状態に応じて定められた部品交換を行う|劣化過程を止めて、修理内容に応じた状態に戻す|
 +|大修理|点検・診断の後、劣化している部品やコンポーネントをすべて新しい部品やコンポーネントと取り替える|劣化過程を止めて、新品同様の状態に戻す|
 +|取替|新品と交換する|新品の状態に戻す|
 +
 +|略語(他動詞)|Keep|Repair|Overhaul|Replace|
 +|状態|稼働を見守る費用|少し修理する費用|新品同様まで修理する費用|新品と交換する費用|
 +|1|C(稼働)[1]|C(小修理)[1]|C(大修理)[1]|C(取替)[1]|
 +|2|C(稼働)[2]|C(小修理)[2]|C(大修理)[2]|C(取替)[2]||
 +|...| | | | |
 +|N|C(稼働)[N]|C(小修理)[N]|C(大修理)[N]|C(取替)[N]|
 +
 +状態を１戻す修理(小修理)のに、2->1と3->2で同じ部品を交換するとも限らない。そのため、修理による状態推移2->1と別の修理3->2同じ費用がかかるとは限らない。
 +
 +また状態を１に戻す修理(大修理)では、2->1と3->1で通常は前者の方が費用が安い。
 +
 +取替の費用は、現在の状態のシステムを下取りに出したり、処分して上がる利益がなければ、一定となる。
 +
 +例えば、次のようなコスト設定があるとする。
 +
 +| |Keep|Repair|Overhaul|Replace|
 +|1|0|10|20|150|
 +|2|0|50|70|150|
 +|3|2000|250|180|150|
 +
 +これをRでは、例えば次のようなコードで与える。
 +
 +<code>
 +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")
 +</code>
 +
 +<code>
 +  Keep  Repair Overhaul Replace
 +1     0     10       20     150
 +2     0     50       70     150
 +3   2000    250      180     150
 +</code>
 +
 +グラフは次のコードで描ける。
 +
 +<code>
 +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"))
 +</code>
 +
 +{{:r:maintenance:costs.png|}}
 +
 +グラフからは、状態が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>
 +
 +