Factorio Space Age 品質上げ方法の評価例（Rによる解析）
この記事はFactorioに関連するものですが、行列計算に関する話が主になりますのでご注意ください。理数系が得意な人で、高校～大学教養課程程度の内容になります。数値計算にはRを用い、使用したコードも掲載しています。
万が一間違いを見つけた場合は教えて下さい。また、よろしければ色んなパターンで計算した結果を教えてくださいね。
前回記事のおさらい
以下の記事の「公式wikiのQuality解説を補足」の続きです。前回の記事を要約すると、品質上げをする際、赤モジュールが使える場合に、赤モジュールと品質モジュールどちらを使うべきかは、自明な解はなく、作成（リサイクル）ルートやモジュールスロット数、prodボーナスで異なるので、個別に評価した方が良い、ということでした。では、特定の例で評価してみましょう。
Wikiの例を検算
まずは、Optimal module usage, Derivationの部分をRで検算します。
定義
計算がしやすいよう、procは汎用な関数にしました。作成、リサイクルどちらでも使用できます。
intv <- matrix( c(1,0,0,0,0),nrow=1,ncol=5)
zerv <- matrix( 0,nrow=1,ncol=5)
zero <- matrix( 0,nrow=5,ncol=5)
exlg <- matrix( c(1,1,1,1,0,1,1,1,1,0,1,1,1,1,0,1,1,1,1,0,1,1,1,1,0),nrow=5,ncol=5)
lege <- matrix( c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1),nrow=5,ncol=5)
proc <- function(p,q) {x <- t( matrix( c((1+p)*(1-q),(1+p)*q*9/10,(1+p)*q*9/100,(1+p)*q*9/1000,(1+p)*q*1/1000,0,(1+p)*(1-q),(1+p)*q*9/10,(1+p)*q*9/100,(1+p)*q*1/100,0,0,(1+p)*(1-q),(1+p)*q*9/10,(1+p)*q*1/10,0,0,0,(1+p)*(1-q),(1+p)*q,0,0,0,0,1+p),nrow=5,ncol=5) ); return(x);}
> intv
[,1] [,2] [,3] [,4] [,5]
[1,] 1 0 0 0 0
> zerv
[,1] [,2] [,3] [,4] [,5]
[1,] 0 0 0 0 0
> zero
[,1] [,2] [,3] [,4] [,5]
[1,] 0 0 0 0 0
[2,] 0 0 0 0 0
[3,] 0 0 0 0 0
[4,] 0 0 0 0 0
[5,] 0 0 0 0 0
> exlg
[,1] [,2] [,3] [,4] [,5]
[1,] 1 1 1 1 1
[2,] 1 1 1 1 1
[3,] 1 1 1 1 1
[4,] 1 1 1 1 1
[5,] 0 0 0 0 0
> lege
[,1] [,2] [,3] [,4] [,5]
[1,] 0 0 0 0 0
[2,] 0 0 0 0 0
[3,] 0 0 0 0 0
[4,] 0 0 0 0 0
[5,] 0 0 0 0 1
> 行列M
以降、Foundary（モジュール数4、prodボーナス+50%）を使用する場合で計算してみます。以下はWikiでの例となったテーブルです。左上がzero、右上が行列P=proc、左下が行列Rはproc*exlg（ただし*はアダマール積）、右下が行列L=legeです。パラメータは品質の確率を正確にしたので若干Wikiとは異なります（Wikiの方が若干不正確です）。
mtrx <- function(pp,pq,rp,rq){ x <- rbind ( cbind (zero,proc(pp,pq)),cbind(proc(rp,rq)*exlg,lege));return(x);}
> mtrx(0.25*0+0.5,0.062*4,-0.75,0.062*4)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0.000 0.0000 0.00000 0.000000 0.0e+00 1.128 0.3348 0.03348 0.003348 0.000372
[2,] 0.000 0.0000 0.00000 0.000000 0.0e+00 0.000 1.1280 0.33480 0.033480 0.003720
[3,] 0.000 0.0000 0.00000 0.000000 0.0e+00 0.000 0.0000 1.12800 0.334800 0.037200
[4,] 0.000 0.0000 0.00000 0.000000 0.0e+00 0.000 0.0000 0.00000 1.128000 0.372000
[5,] 0.000 0.0000 0.00000 0.000000 0.0e+00 0.000 0.0000 0.00000 0.000000 1.500000
[6,] 0.188 0.0558 0.00558 0.000558 6.2e-05 0.000 0.0000 0.00000 0.000000 0.000000
[7,] 0.000 0.1880 0.05580 0.005580 6.2e-04 0.000 0.0000 0.00000 0.000000 0.000000
[8,] 0.000 0.0000 0.18800 0.055800 6.2e-03 0.000 0.0000 0.00000 0.000000 0.000000
[9,] 0.000 0.0000 0.00000 0.188000 6.2e-02 0.000 0.0000 0.00000 0.000000 0.000000
[10,] 0.000 0.0000 0.00000 0.000000 0.0e+00 0.000 0.0000 0.00000 0.000000 1.000000convは、収束したと判定するまで繰り返し行列を掛けて確率計算を行います。これを用いて最適なp,qの組み合わせを調査します。出力内容はループ回数（往復でn+2）およびアイテム数（Nomarl,Legend）です。実際はモジュール数は整数ですが、どこにピークがあるかを探るため整数でない数字でも確認しました。その結果、ちょうど赤モジュール3.0個、品質モジュール1.0個あたりにピークがあることがわかります。また、赤3、品質1は、赤0、品質4に比べて5割Legendの数が大きいです。結構差が出ますね。WikiではFoundaryの場合は25の材料で1つのLegendaryが見込めるとありますので、結果は一致します。
ちなみに、2番目の引数を0にするとエラーになります。記事初掲載後、vctn[9]/vctn[10]の箇所でvctn[10]=0になるためと指摘いただきました。しかし、ここを直すと全体に影響があって大変なので、式はこのまま修正しないこととします。0でない十分小さい数字を使用するなど、運用でカバーしてください。
conv <- function(pp,pq){n <- 1; r <- T; x <- diag(10); while(r){ n <- n+1; x <- x %*% mtrx(pp,pq,-0.75,0.062*4);vctn <- c(intv,zerv) %*% x; r <- (n %% 2) | vctn[9]/vctn[10] > 0.01;};print(n);return(c(vctn[6],vctn[10])); }
> conv(0.5+0.25*0,0.062*4)
[1] 18
[1] 4.613650e-06 2.450029e-02
> conv(0.5+0.25*1,0.062*3)
[1] 20
[1] 1.009358e-05 3.083419e-02
> conv(0.5+0.25*2,0.062*2)
[1] 22
[1] 2.633017e-05 3.582863e-02
> conv(0.5+0.25*2.9,0.062*1.1)
[1] 26
[1] 2.549057e-05 3.799474e-02
> conv(0.5+0.25*3.0,0.062*1.0)
[1] 28
[1] 1.274815e-05 3.807804e-02
> conv(0.5+0.25*3.1,0.062*0.9)
[1] 28
[1] 1.631878e-05 3.804108e-02
> conv(0.5+0.25*4,0.062*0.001)
[1] 34
[1] 0.0000141595 0.0349740478
>
> conv(0.5+0.25*4,0.062*0)
while (r) { でエラー: TRUE/FALSE が必要なところが欠損値です
> Holmium plateの品質上げ
さて、今度はHolmium plateの品質上げについて考えてみます。
名称について
私のゲームは日本語に翻訳していないので、正確な日本語表記がわかりません。この記事では表記は英語の名称を使います。英語の正式名称がわかるよう念のためキャプチャを貼っておきます。
問題設定
Holmium plateは、リサイクラーにかけると自分自身に戻ります。リサイクルを繰り返して品質上げを行うのはとても資源効率が悪いです。リサイクルのみの繰り返しは、その資源の価値が低く十分高速に供給できるか、他に方法がない場合のみ行うべきです。Holmium oreは希少な資源です。
そこで、以下のルートでHolmium plateの品質上げをする場合を考えます。赤モジュールと品質モジュールの選択の余地があるのは(*1)だけです。また、EleMag plantを作成する際、Holmium plate以外の材料は十分に使用可能であるとします。
Holmium solution
→Foundary(Holmium solution)=Holmium plate (*1)
→EleMag plant(Holmium plate)=EleMag plant
→Recycler(EleMag plant)=Holmium plate
→EleMag plant(Holmium plate)=EleMag plant
→Recycler(EleMag plant)=Holmium plate
︙
→EleMag plant(Holmium plate)=EleMag plant
→Recycler(EleMag plant)=Holmium plate
式はこんな感じになります。リサイクルが最終工程になります。式の共通化や理解をしやすくするため、Wikiの例と比べて、各ベクトルの要素は左右交換、行列Mは右上と左下を交換します。
Rによる解析
M_f=mtrf、M_em=mtrxです。mtrxは前回と比べて引数の指定順が異なりますが、mtrx自体は同じです。なお、mtrfは6行目しか計算に使用されません。
mtrf <- function(pp,pq){ x <- rbind ( cbind (zero,zero),cbind(zero,proc(pp,pq)));return(x);}
mtrx <- function(pp,pq,rp,rq){ x <- rbind ( cbind (zero,proc(pp,pq)),cbind(proc(rp,rq)*exlg,lege));return(x);}
> mtrf(0.5,0.062*4)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0 0 0 0 0 0.000 0.0000 0.00000 0.000000 0.000000
[2,] 0 0 0 0 0 0.000 0.0000 0.00000 0.000000 0.000000
[3,] 0 0 0 0 0 0.000 0.0000 0.00000 0.000000 0.000000
[4,] 0 0 0 0 0 0.000 0.0000 0.00000 0.000000 0.000000
[5,] 0 0 0 0 0 0.000 0.0000 0.00000 0.000000 0.000000
[6,] 0 0 0 0 0 1.128 0.3348 0.03348 0.003348 0.000372
[7,] 0 0 0 0 0 0.000 1.1280 0.33480 0.033480 0.003720
[8,] 0 0 0 0 0 0.000 0.0000 1.12800 0.334800 0.037200
[9,] 0 0 0 0 0 0.000 0.0000 0.00000 1.128000 0.372000
[10,] 0 0 0 0 0 0.000 0.0000 0.00000 0.000000 1.500000
> convは前回と少し違います。計算すると、論理的には赤モジュール2.5個、品質モジュール1.5個のあたりがピークで、実際のゲームでは赤3質1が一番資源効率が良いことがわかります。ただ、赤0質4の最も悪いケースと比べて1割しか違いはありません。どのケースでもそれほど差はないように感じます。
conv <- function(pp,pq){n <- 1; r <- T; x <- mtrf(pp,pq); while(r){ n <- n+1; x <- x %*% mtrx(-0.75,0.062*4,0.5+0.25*0,0.062*5);vctn <- c(zerv,intv) %*% x; r <- !(n %% 2) | vctn[9]/vctn[10] > 0.01;};print(n);return(c(vctn[6],vctn[10])); }
> conv(0.5+0.25*0,0.062*4)
[1] 15
[1] 1.191232e-05 3.040103e-02
> conv(0.5+0.25*1,0.062*3)
[1] 15
[1] 1.504353e-05 3.225977e-02
> conv(0.5+0.25*2,0.062*2)
[1] 17
[1] 3.600143e-06 3.327743e-02
> conv(0.5+0.25*2.4,0.062*1.6)
[1] 17
[1] 3.887168e-06 3.340229e-02
> conv(0.5+0.25*2.5,0.062*1.5)
[1] 17
[1] 3.960516e-06 3.341060e-02
> conv(0.5+0.25*2.6,0.062*1.4)
[1] 17
[1] 4.034502e-06 3.340975e-02
> conv(0.5+0.25*3,0.062*1)
[1] 17
[1] 4.336816e-06 3.331475e-02
> conv(0.5+0.25*4,0.062*0)
[1] 17
[1] 5.13719e-06 3.24360e-02
>