差分

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

この比較画面へのリンク

両方とも前のリビジョン前のリビジョン
次のリビジョン		前のリビジョン
r:maintenance:condition_monitoring [2018/12/17 00:44] – [費用関数の設定] 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>
  
-==== 遷移行列の配列の設定 ==== 
  
-MDPtoolboxパッケージのために、行動ごとの遷移行列をすべて１つの配列に収める。 
-次のようにオブジェクトPを設定する。 
- 
-<code> 
-P = array(dim=c(n.status, n.status, n.action)) 
-P[,,1] = P.Dgr 
-P[,,2] = P.Rpr 
-P[,,3] = P.Rpl 
-</code> 
- 
- 
-==== 費用関数の設定 ==== 
  
-同じく、MDPtoolboxパッケージのために費用関数を設定する。 +==== MDPtoobox利用の準備 ====
-行動の順序は、上の遷移行列の配列と同じにする必要がある。+
  
 マルコフ決定過程の定義に必要なものは次の２つ。 マルコフ決定過程の定義に必要なものは次の２つ。
行 559: 行 578:
   * 行動と状態ごとの費用行列R   * 行動と状態ごとの費用行列R
  
-Rは手作業での設定になるので、ついでにPもこのページの上の方で用いた状態が3種類、行動が3種類ある場合の例を設定してみる。この遷移行列は3×3×3の3次元配列となる。+Rは手作業での設定になる。
  
-<code> +===== 遷移行列の配列の設定 =====
-n.status = 3 +
-n.action = 3 +
-P = array(dim=c(n.status,n.status,n.action)) +
-</code>+
  
-配列Pに行動ごとの遷移行列を代入していく。 +MDPtoolboxパッケージのために、行動ごとの遷移行列をすべて１つの配列に収める。 
- +次のようにオブジェクトPを設定する。
-<code> +
-P[,,1] = matrix(c( +
-  9/10,1/10,0, +
-  0,9/10,1/10, +
-  0,0,1),nrow=3,ncol=3,byrow=TRUE) +
-P[,,2] = matrix(c( +
-  1,0,0, +
-  1,0,0, +
-  0,1,0),nrow=3,ncol=3,byrow=TRUE) +
-P[,,3] = matrix(c( +
-  1,0,0, +
-  1,0,0, +
-  1,0,0),nrow=3,ncol=3,byrow=TRUE) +
-</code> +
- +
-上の方に記したコードを実行してあれば、次のように代入してもいい。+
  
 <code> <code>
 +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
行 592: 行 593:
 </code> </code>
  
-できあがった配列Pの中身を確認してみる。+P.Dgrを正しく推定できていれば、Pはn.status×n.status×n.actionの3次元配列になる。
  
-<code> +===== 費用関数の設定 =====
-P +
-</code>+
  
-次に費用関数を定義する。+同じく、MDPtoolboxパッケージのために費用関数を設定する。ここは、以下のコードをそのまま貼るのではなく、数字の数を状態数に合わせて適宜変更する必要がある。行動の数n.actionは3としている。また行動の順序は、上の遷移行列の配列と同じにする必要がある。
  
 <code> <code>
-R = matrix(c( +C.Dgr = c(0,0,2000)
-    0,10,150, +
-    0,50,150, +
-    2000,250,150),nrow=n.status, ncol=n.action, byrow=TRUE) +
-colnames(R) = c("keep","repair","replace") +
-rownames(R) = c(1:3) +
-R = -R +
-</code> +
- +
-これも先ほどのコードを実行してあれば、次のようにしてもいい。 +
-<code> +
-C.Opr = c(0,0,2000)+
 C.Rpr = c(10,50,250)  C.Rpr = c(10,50,250) 
 C.Rpl = c(150,150,150) C.Rpl = c(150,150,150)
-Cost = cbind(C.Opr, C.Rpr, C.Rpl)+Cost = cbind(C.Dgr, C.Rpr, C.Rpl)
 colnames(Cost) = c("Keep","Repair","Replace") colnames(Cost) = c("Keep","Repair","Replace")
 rownames(Cost) = c("1","2","3") rownames(Cost) = c("1","2","3")
行 621: 行 609:
 </code> </code>
  
-できあがった配列Pの中身を確認してみる。+6状態であれば次のようにする。
  
 <code> <code>
-R+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>
  
 +各自の費用関数は、下記の数値の数を入れ替えて設定することになる。
 +
 +<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に設定して、最適方策を試算している。
 === 価値反復法 === === 価値反復法 ===