| 図８．１．１ 机上実験 |
 |
| | HOME | 目次 | | 2009 (C) H.Ishikawa |
| ８．１ | 机上で実験してみる | ||
| ８．２ | コンピュータによる例題のシミュレーション－－固定時間方式－－ | ||
| ８．３ | 高速な方法－－変化点方式－－ | ||
| ８．４ | どちらのシミュレーション方式がいいか | ||
| ８．５ | シミュレーション時間の決め方 | ||
| 演習 |
前の章までで、でいくつかの簡単なシミュレーションを試みましたが、それらはコインを投げるとか、針を落とすとか、同時には１つの事象しか起こらない場合でした。しかし、一般には同時に並行的に、いくつもの事象が進行する場合がたくさんあります。この同時進行のシミュレーションを行なうにはどうしたらいいでしょうか。この章では、機械修理の問題という、複数の機械があり、複数の修理者がいるという、同時進行の典型的なものをとりあげ、同時進行のシミュレーション技法をマスタすることとします。
８．１ 机上で実験してみる
機械修理問題とは
今ここに、m 台の生産機械があるとします。経験からすると、平均 a 時間すると具合が悪くなり故障します。その修理調整は１人の修理者によって
h 時間かかります。生産機械は１時間あたり r 円の収益（収入）があり、修理者は１人あたり１時間
c 円の人件費がかかるものとします。何人の修理者を用意し修理させると、利益（収益マイナス費用）が最大となるか、これが機械修理の問題です1)。
故障発生間隔、修理時間それぞれが指数分布の場合、理論的解析が可能ですが2)、ここではこれをシミュレーションにより解析してみましょう。
注
1) スウェーデンの工業で応用され成功しているものとして、文献(1)に紹介がある。
2) ８．４節に理論解析のプログラムを示す。
机上実験
この例題のシミュレーションは機械 m 台、修理者n人が、独立な動きをしますから、そう簡単ではありません。コンピュータを用いる前に、机上実験をカードを用いて試みて同時動作の概念をつかむこととしましょう。
機械の台数 m=5 台、修理者の人数 n=2 人、故障の平均発生間隔 a=10 時間、平均の修理時間 h=4 時間と仮定します。
＜準備するもの＞
故障発生機構：１０枚のカード。そのうちの１枚に「故障発生」と書いておく（１０時間に１回故障することを、意味する）。
修理完了機構：４枚のカード。そのうちの１枚に「修理完了」と書いておく（修理時間が４時間を意味する）。
機械カード：５枚のカード。機械の番号１～５を書いておく。
記録用紙：シミュレーションの経過を記録する用紙。
＜机上シミュレーション＞
シミュレーションの手順は次のとおりです。
① 初期設定
図８．１．１のように机の上の「機械」の位置に機械カードを５枚ならべる。また「待ち行列」の位置、「修理者」の位置を空けておく。記録紙に、時刻、機械の状態、修理者が何番の機械の保守をしているかを記入できるようにしておく。時刻は第１時間目と書いておく。
② 故障発生
機械１のカードが「機械」の位置においてあれば、それは稼働中を意味します。「故障発生機構」のカードをよくきって、そのうちの１枚を引く。「故障発生」のカードを引いたら、機械１のカードを「待ち行列」の位置にならべる。記録用紙の１時間目の機械１の欄に×（故障発生のマーク）をつける。
③ 故障発生繰返し
機械２から５について②を繰り返す。
④ 修理開始
次に修理者１の位置に機械カードがなければ、その修理者は修理中でない未稼働の意味であるので、「待ち行列」の先頭の機械カードを取り出し、修理者１のところにおく。そして、記録用紙の修理者１の欄に、修理する機械番号を１と書いておく。
⑤ 修理完了
もしも修理者１の位置に、機械カードが置いてあれば、それは修理者が修理中であること、つまり稼働中があることを意味している。「修理完了機構」のカードをよくきってその中から１枚を引く。もし「修理完了」カードを引いたら、その時点で修理完了であるから、「修理者」の位置にあった機械カードをもとの「機械」の位置にもどす。記録用紙の機械の欄に○（修理完了のマーク）をつける。
⑥ 修理完了繰返し
修理者２について④⑤を繰り返す。
⑦ 以上の事が、１時間の中で同時に起った事です。時計を１時間進め、②から繰り返す。
以上により、数十時間分を繰り返し、記録用紙を見れば、表８．１．１のようになるでしょう。各機械が何時間故障であったか、修理開始まで何時間待たされるか、修理者はどのくらい稼働していたか、この表から計算することができるでしょう。
この方法は、３．１節で述べたベルヌーイ試行の実験の応用で、時計を一時停めて、すべての機械、修理者に対し、ベルヌーイ試行を行った後、時計を進めると考えて下さい。これが同時進行のシミュレーション方法です。
| 図８．１．１ 机上実験 |
|
| 表８．１．１ 記録用紙 |
 |
| この章始め |
| 図８．２．１ 機械修理問題の状態遷移図 |
 |
| 図８．２．２ 待ち行列の構造 |
 |
| 図８．２．３ 固定時間方式シミュレーション結果 (部分） |
 |
| S0820.java 機械修理の問題 (固定時間方式) | Download | 実行 | |
/* S0820.java
* 機械修理の問題 (固定時間方式)
* (C) H.Ishikawa 2008
*/
package simulation;
import java.applet.*;
import java.awt.*;
import java.text.DecimalFormat;
public class S0820 extends Applet{
int M = 5; /* 機械の台数 */
int N = 2; /* 修理者の人数 */
double A = 10.0; /* 平均故障発生間隔 */
double H = 4.0; /* 平均修理時間 */
double DT = 0.5; /* 時刻増分 */
double T_END = 30.0; /* シミュレーションを終る時刻 */
public void paint(Graphics g) {
double tt = 0.0; /* シミュレーションクロック */
int i; /* 待ち行列のforのカウンタ */
int j; /* 機械のforのカウンタ */
int k; /* 修理者のforのカウンタ */
int f = 0; /* 故障している機械の台数 */
int q = 0; /* 待ち行列の長さ */
int l = 0;
int sm[] = new int[M + 1]; /* 機械jの状態 */
int sn[] = new int[N + 1]; /* 修理者kの状態 */
int sq[] = new int[M + 1]; /* 待ち行列 */
double p[] = new double[M + 1]; /* 機械がf台故障している時間の割合 */
DecimalFormat exFormat1 = new DecimalFormat(" 0.000000");
g.drawString(" tt sm[j] sn[k] sq[i]", 10, 10);
/* 開始 */
while (tt <= T_END) {
/* (1)故障発生 */
for (j = 1; j <= M; j ++) {
if (sm[j] == 0) { /* 機械が稼働中なら */
if ((1.0 / A * DT) > Math.random()) { //(2)
sm[j] = 1; /* ある確率で故障する */
q = q + 1; /* (3)待ち行列の後ろにつく */
sq[q] = j;
}
}
}
/* (4)故障修理 */
for (k = 1; k <= N; k ++) {
if (sn[k] == 0) { /* (5)修理者が未稼働なら */
if (q > 0) {
sn[k] = sq[1]; /* 待ち行列から取り出す */
q = q - 1;
if (q > 0) {
for (i = 1; i <= q; i ++) {
sq[i] = sq[i + 1]; /* (6)待ち行列を1進める */
}
}
}
} else {
if ((1.0 / H * DT) > Math.random()) { /* (7)修理者は機械を修理し */
sm[sn[k]] = 0; /* 稼働状態にする */
sn[k] = 0;
}
}
}
/* (8)表示 */
g.drawString(" " + tt, 10, 20 + 10 * l);
for (j = 1; j <= M; j ++) g.drawString(" " + sm[j], 50 + 10 * j, 20 + 10 * l );
for (k = 1; k <= N; k ++) g.drawString(" " + sn[k], 100 + 10 * M +10 * k, 20 + 10 * l );
for (i = 1; i <= q; i ++) g.drawString(" " + sq[i], 150+ 10 * (M + N) + 10 * i, 20 + 10 * l );
l= l + 1;
/* 故障している機械の台数を数える */
f = 0;
for (j = 1; j <= M; j ++) {
if (sm[j] == 1) { f = f + 1; }
}
p[f] = p[f] + DT;
/* (9)時刻を進める */
tt =tt + DT;
}
/* 結果出力 */
for (j = 0; j <= M; j ++) {
g.drawString("p["+ j + "]="+ " " + exFormat1.format((double)p[j] / T_END) ,10 , 20 + 10 * l + 10 * j);
}
}
}
|
| この章始め |
| 図８．３．１ 固定時間方式と変化点方式 |
 |
| 図８．３．２ 変化点方式シミュレーション結果 |
 |
| S0830.java 機械修理の問題 (変化点方式) | Download | 実行 | |
/* S0830.java
* 機械修理の問題 (変化点方式)
* (C) H.Ishikawa 2008
*/
package simulation;
import java.applet.*;
import java.awt.*;
import java.text.DecimalFormat;
public class S0830 extends Applet{
int M = 5; /* 機械の台数 */
int N = 2; /* 修理者の人数 */
double A = 10.0; /* 平均故障発生間隔 */
double H = 4.0; /* 平均修理時間 */
double DT = 0.5; /* 時刻増分 */
double T_END = 100.0; /* シミュレーションを終る時刻 */
double INF = 1.0e100; /* 十分大きい値 */
double EPS = 1.0e-100; /* 十分小さい値 */
public void paint(Graphics g) {
double tt = 0.0; /* シミュレーションクロック */
int i; /* 待ち行列のforのカウンタ */
int j; /* 機械のforのカウンタ */
int k; /* 修理者のforのカウンタ */
int f = 0; /* 故障している機械の台数 */
int q = 0; /* 待ち行列の長さ */
int l = 0; /* 表示行 */
int sm[] = new int[M + 1]; /* 機械jの状態 */
double tm[] = new double[M + 1]; /* (1)機械に次の状態変化が起るまでの時間 */
int sn[] = new int[N + 1]; /* 修理者kの状態 */
double tn[] = new double[N + 1]; /* (1)修理者に次の状態変化が起るまでの時間 */
int sq[] = new int[M + 1]; /* 待ち行列 */
double p[] = new double[M + 1]; /* (1)機械がf台故障している時間の割合 */
double at; /* 時刻増分 */
DecimalFormat exFormat1 = new DecimalFormat(" 0.000000");
DecimalFormat exFormat2 = new DecimalFormat(" 0.00");
g.drawString(" tt sm[j] sn[k] sq[i]", 10, 10);
/* (2)初期設定 */
for (j = 1; j <= M; j ++) {
sm[j] = 0; tm[j] = -A * Math.log(1 - Math.random());
}
for (k = 1; k <= N; k ++) {
sn[k] = 0; tn[k] = INF;
}
/* (3)開始 */
while (tt <= T_END) {
/* 故障発生 */
for (j = 1; j <= M; j ++) {
if (Math.abs(tm[j]) < EPS) { /* (4)機械jが故障発生 */
sm[j] = 1; tm[j] = INF; /* (5) */
q = q + 1;
sq[q] = j;
}
}
/* (6)故障修理 */
for (k = 1; k <= N; k ++) {
if (Math.abs(tn[k]) < EPS) { /* (7)修理者kの修理が完了 */
sm[sn[k]] = 0; tm[sn[k]] = -A * Math.log(1 - Math.random()); /* (8) */
}
if (Math.abs(tn[k]) <EPS || sn[k] == 0) { /* (9)kが未稼動の場合 */
if (q > 0) { /* (10)待ち行列を取り出す */
sn[k] = sq[1]; tn[k] = -H * Math.log(1 - Math.random());
q = q - 1;
if (q > 0) {
for (i = 1; i <= q; i ++) {
sq[i] = sq[i + 1];
}
}
} else { /* 待ち行列が空なら修理者kは未稼働状態になる */
sn[k] = 0; tn[k] = INF;
}
}
}
/* (13)表示 */
g.drawString(" " + exFormat2.format(tt), 10, 20 + 10 * l);
for (j = 1; j <= M; j ++) g.drawString(" " + sm[j], 50 + 10 * j, 20 + 10 * l );
for (k = 1; k <= N; k ++) g.drawString(" " + sn[k], 100 + 10 * M +10 * k, 20 + 10 * l );
for (i = 1; i <= q; i ++) g.drawString(" " + sq[i], 150+ 10 * (M + N) + 10 * i, 20 + 10 * l );
l= l + 1;
/* (11)次に状態変化の起る時刻を求める */
at = INF;
for (j = 1; j <= M; j ++) {
if (tm[j] < at) at = tm[j];
}
for (k = 1; k <= N; k++) {
if (tn[k] < at) at = tn[k];
}
/* (14)故障している機械の台数を数える */
f = 0;
for (j = 1; j <= M; j ++) {
if (sm[j] == 1) { f ++; }
}
p[f] = p[f] + at;
/* (12)時間を進める */
for (j = 1; j <= M; j ++) {
tm[j] = tm[j] - at;
}
for (k = 1; k <= N; k ++) {
tn[k] = tn[k] - at;
}
tt = tt+ at;
}
/* 結果出力 */
for (j = 0; j <= M; j ++) {
g.drawString("p["+ j + "]="+ " " + exFormat1.format((double)p[j] / T_END) ,10 , 20 + 10 * l + 10 * j);
}
}
}
|
| この章始め |
８．４ どちらのシミュレーション方式がいいか
| 図８．４．１ |
 |
 |
・・・ | 式８．４．１ |
 |
・・・・・ | 式８．４．２式 |
 |
・・・ | 式８．４．３ |
 |
・・・ | 式８．４．４ |
| 図８．４．２ S0840.javaの 結果 |
 |
| 表８．４．１ 固定時間方式と変化点方式の比較 |
 |
| S0840.java 機械修理の問題 (理論値) | Download | 実行 | |
/* S0840.java
* 機械修理の問題 (理論値)
* (C) H.Ishikawa 2008
*/
package simulation;
import java.applet.*;
import java.awt.*;
import java.text.DecimalFormat;
public class S0840 extends Applet{
int M = 5; /* 機械の台数 */
int N = 2; /* 修理者の人数 */
double A = 10.0; /* 平均故障発生間隔 */
double H = 4.0; /* 平均修理時間 */
public void paint(Graphics g) {
int j; /* 機械のforのカウンタ */
double p[] = new double[M + 1]; /* 機械がf台故障している時間の割合 */
double p0 = 0.0;
DecimalFormat exFormat1 = new DecimalFormat(" 0.000000");
p[0] = 1;
p[1] = (double)M * H / A * p[0];
for (j = 1; j <= N - 1; j ++) {
p[j + 1] = (double)(M - j) / (j + 1) * H / A * p[j];
}
for (j = N; j <= M - 1; j ++) {
p[j + 1] = (double)(M - j) / N * H / A * p[j];
}
for (j = 0; j<= M; j ++) {
p0 = p0 + p[j];
}
for (j = 0; j<= M; j ++) {
p[j] = p[j] / p0;
}
for (j = 0; j <= M; j ++) {
g.drawString("p["+ j + "] ="+ " " + exFormat1.format(p[j] ) ,10 , 20 + 10 * j);
}
}
}
|
| この章始め |
８．５ シミュレーション時間の決め方
シミュレーションをおこなう場合，信頼に足る結果を得るのにどのくらい計算機を回すのが必要でしょうか。JAVAアプレットを使えば、最終結果をアウトプットするのみならず，途中の経過をディスプレイすることができますから，その形をみて，収束状況を，判定することができます。
ディスプレイの画面には，横軸に試行回数あるいはシミュレーションクロック，縦軸に結果をとると，図８．５．１のようなグラフが表示されるでしょう。すなわち，始めはばらつきが大きく，しだいに，収束されて行きます。３．４節の大数の法則と，３．５節の中心極限定理を，思いだして下さい。縦軸はある試行の結果，成功の回数と試行の回数の比ですから，大数の法則により，試行回数
n を多くすると一定値に近づくわけです。しからば，この場合のバラツキはどうでしょうか。中心極限定理は， x の分布にかかわらず， x を n 個合計し， n で割ると，標準偏差が， １／√n の正規分布となるといっています。 n を大きくすると標準偏差は小さく中心に近ずくわけです。しかし，その割合は n を，１００倍しても，標準偏差は，１／１０にしかなりません。つまり，１桁精度を上げようとすれば，シミュレーション時間を１００倍かける必要があります。
| 図８．５．１ シミュレーション時間と結果の収束 |
 |
| この章始め |
１． 機械の台数＝５台、故障の平均発生間隔＝１０時間、平均修理時間＝４時間、１時間あたりの機械１台の収益＝８万円、修理者１人１時間あたりの費用＝２万円としたとき、最大利益が得られる理論的な修理者の人数を求めなさい。
（解答） A0810.java | 実行 |
これによれば、２人が最良
２． プログラムS0830.javaで、修理時間を位相k=3のアーラン分布とした場合について、シミュレーションしなさい。
（解答） A0820.java | 実行 |
３． プログラムS0830.javaを改造して、故障している機械の割合 pf[j]/tt を、時々刻々表示し、シミュレーションの収束状況を、観測できるようにしなさい。また機械の稼働率、修理者の稼働率、機械の待ち時間の分布、機械の待ち合わせ平均時間、その分散を求めるようにしなさい。
（解答） A0830.java | 実行 |
４． １本の動力線を m 人の溶接工が共用しています。溶接工はランダムに断続的に１Ａの電流を使用します。その時間の平均値は h 分で間隔の平均は a 分です。電流の波形をシミュレーションしなさい。
（解答） A0840.java | 実行 |
| この章始め |
| | HOME | 目次 | | 2009 (C) H.Ishikawa |