差分

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

この比較画面へのリンク

両方とも前のリビジョン前のリビジョン
次のリビジョン		前のリビジョン
r:maintenance:condition_monitoring [2018/12/16 18:00] – 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>
 +
 +
 +=== 劣化と保全 ===
  
 故障時間のデータしかなけば、寿命分布を推定し、時間取り替えやブロック取り替えを考えるしかない。 故障時間のデータしかなけば、寿命分布を推定し、時間取り替えやブロック取り替えを考えるしかない。
行 383: 行 401:
 plot(c(1,n.epoch),c(min(Data),max(Data)), plot(c(1,n.epoch),c(min(Data),max(Data)),
      xlab="Time",ylab="Data",type="n")      xlab="Time",ylab="Data",type="n")
-for( i in c(1:n.b) ) {+for( i in c(1:n.item) ) {
   lines(c(1:n.epoch),Data[i,])   lines(c(1:n.epoch),Data[i,])
 } }
 +</code>
 +
 +別のグラフの描き方もある。
 +
 +<code>
 +matplot(t(Data))
 </code> </code>
  
行 395: 行 419:
 </code> </code>
  
-マルコフ決定過程で保全方策を最適化するには、このグラフのデータを離散の状態の間の推移に変換する必要がある。等間隔の区間で離散化するために、最小値status.min、最大値status.max、区間幅status.widthを定める。+マルコフ決定過程で保全方策を最適化するには、このグラフのデータを離散の状態の間の推移に変換する必要がある。等間隔の区間で離散化するために、最小値status.min、最大値status.max、区間幅status.widthを定める。次の数値は、例であり、各自で設定する必要がある。
  
 <code> <code>
行 401: 行 425:
 status.max = 1.05 status.max = 1.05
 status.width = 0.1 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> </code>
  
行 406: 行 438:
  
 <code> <code>
-n.status = (status.max-status.min)/tic.diff+n.status = ceiling((status.max-status.min)/status.width)
 </code> </code>
  
行 431: 行 463:
 plot(c(1,n.epoch),c(min(Data.states),max(Data.states)), plot(c(1,n.epoch),c(min(Data.states),max(Data.states)),
      xlab="Time",ylab="Status",type="n")      xlab="Time",ylab="Status",type="n")
-for( i in c(1:n.b) ) {+for( i in c(1:n.item) ) {
   lines(c(1:n.epoch),Data.states[i,])   lines(c(1:n.epoch),Data.states[i,])
 } }
行 463: 行 495:
 for( i in c(1:n.status) ) { for( i in c(1:n.status) ) {
   if(sum((Data>status.breaks[i]) & (Data<=status.breaks[i+1]))>0) {   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)+    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)), plot(c(1,n.epoch),c(min(Data.states),max(Data.states)),
      xlab="Time",ylab="Status",type="n")      xlab="Time",ylab="Status",type="n")
-for( i in c(1:n.b) ) {+for( i in c(1:n.item) ) {
   lines(c(1:n.epoch),Data.states[i,])   lines(c(1:n.epoch),Data.states[i,])
 } }
行 510: 行 542:
  
 <code> <code>
-rbind(c(1,rep(0,n.status-1)), +n.translation = n.status-1 
-            cbind(diag(1,nrow=n.status-1,ncol=n.status-1),0))+rbind(c(1,rep(0,n.status)), 
 +            cbind(diag(rep(1,n.status)),0))
 </code> </code>
  
行 517: 行 550:
 <code> <code>
 P.Rpr = rbind(c(1,rep(0,n.status-1)), P.Rpr = rbind(c(1,rep(0,n.status-1)),
-            cbind(diag(1,nrow=n.status-1,ncol=n.status-1),0))+            cbind(diag(rep(1,n.status-1)),0))
 rownames(P.Rpr) = c(1:n.status) rownames(P.Rpr) = c(1:n.status)
 colnames(P.Rpr) = c(1:n.status) colnames(P.Rpr) = c(1:n.status)
行 525: 行 558:
  
 <code> <code>
-matrix(rep(c(1,rep(0,n.status-1)),n.status),nrow=n.status,ncol=n.status,byrow=TRUE)+cbind(1,matrix(0,n.status,n.status-1))
 </code> </code>
  
行 531: 行 564:
  
 <code> <code>
-P.Rpl = matrix(rep(c(1,rep(0,n.status-1)),n.status),nrow=n.status,ncol=n.status,byrow=TRUE)+P.Rpl = cbind(1,matrix(0,n.status,n.status-1))
 rownames(P.Rpl) = c(1:n.status) rownames(P.Rpl) = c(1:n.status)
 colnames(P.Rpl) = c(1:n.status) colnames(P.Rpl) = c(1:n.status)
 </code> </code>
  
-==== 遷移行列の配列の設定 ====+ 
 + 
 +==== MDPtoobox利用の準備 ==== 
 + 
 +マルコフ決定過程の定義に必要なものは次の２つ。 
 + 
 +  * 行動ごとの遷移行列(推移確率行列)P 
 +  * 行動と状態ごとの費用行列R 
 + 
 +Rは手作業での設定になる。 
 + 
 +===== 遷移行列の配列の設定 =====
  
 MDPtoolboxパッケージのために、行動ごとの遷移行列をすべて１つの配列に収める。 MDPtoolboxパッケージのために、行動ごとの遷移行列をすべて１つの配列に収める。
行 542: 行 586:
  
 <code> <code>
-P = array(dim=c(n.status, n.status, n.action)+n.action = 3 
 +P = array(dim=c(n.status, n.status, n.action))
 P[,,1] = P.Dgr P[,,1] = P.Dgr
 P[,,2] = P.Rpr P[,,2] = P.Rpr
行 548: 行 593:
 </code> </code>
  
 +P.Dgrを正しく推定できていれば、Pはn.status×n.status×n.actionの3次元配列になる。
  
-==== 費用関数の設定 ====+===== 費用関数の設定 =====
  
-同じく、MDPtoolboxパッケージのために費用関数を設定する。 +同じく、MDPtoolboxパッケージのために費用関数を設定する。ここは、以下のコードをそのまま貼るのではなく、数字の数を状態数に合わせて適宜変更する必要がある。行動の数n.actionは3としている。また行動の順序は、上の遷移行列の配列と同じにする必要がある。
-行動の順序は、上の遷移行列の配列と同じにする必要がある。+
  
 <code> <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> </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>
 +