 |
|
８．同時進行のシミュレーション技法－－機械修理の問題を例題に－－
|
前の章までで、でいくつかの簡単なシミュレーションを試みましたが、それらはコインを投げるとか、針を落とすとか、同時には１つの事象しか起こらない場合
でした。しかし、一般には同時に並行的に、いくつもの事象が進行する場合がたくさんあります。この同時進行のシミュレーションを行なうにはどうしたらいい
でしょうか。この章では、機械修理の問題という、複数の機械があり、複数の修理者がいるという、同時進行の典型的なものをとりあげ、同時進行のシミュレー
ション技法をマスタすることとします。
機械修理問題とは
今ここに、m 台の生産機械があるとします。経験からすると、平均 a 時間すると具合が悪くなり故障します。その修理調整は１人の修理者によって
h 時間かかります。生産機械は１時間あたり r 円の収益（収入）があり、修理者は１人あたり１時間
c 円の人件費がかかるものとします。何人の修理者を用意し修理させると、利益（収益マイナス費用）が最大となるか、これが機械修理の問題です1)。
故障発生間隔、修理時間それぞれが指数分布の場合、理論的解析が可能ですが2)、ここではこれをシミュレーションにより解析してみましょう。
注
1) スウェーデンの工業で応用され成功しているものとして、文献(1)に紹介がある。
2) ８．４節に理論解析のプログラムを示す。
机上実験
この例題のシミュレーションは機械 m 台、修理者n人が、独立な動きをしますから、そう簡単ではありません。コンピュータを用いる前に、机上実験をカードを用いて試みて同時動作の概念をつかむこととしましょう。
機械の台数 m=5 台、修理者の人数 n=2 人、故障の平均発生間隔 a=10 時間、平均の修理時間 h=4 時間と仮定します。
＜準備するもの＞
故障発生機構：１０枚のカード。そのうちの１枚に「故障発生」と書いておく（１０時間に１回故障することを、意味する）。
修理完了機構：４枚のカード。そのうちの１枚に「修理完了」と書いておく（修理時間が４時間を意味する）。
機械カード：５枚のカード。機械の番号１～５を書いておく。
記録用紙：シミュレーションの経過を記録する用紙。
＜机上シミュレーション＞
シミュレーションの手順は次のとおりです。
① 初期設定
図８．１．１のように机の上の「機械」の位置に機械カードを５枚ならべる。また「待ち行列」の位置、「修理者」の位置を空けておく。記録紙に、時刻、機械の状態、修理者が何番の機械の保守をしているかを記入できるようにしておく。時刻は第１時間目と書いておく。
② 故障発生
機械１のカードが「機械」の位置においてあれば、それは稼働中を意味します。「故障発生機構」のカードをよくきって、そのうちの１枚を引く。「故障発
生」のカードを引いたら、機械１のカードを「待ち行列」の位置にならべる。記録用紙の１時間目の機械１の欄に×（故障発生のマーク）をつける。
③ 故障発生繰返し
機械２から５について②を繰り返す。
④ 修理開始
次に修理者１の位置に機械カードがなければ、その修理者は修理中でない未稼働の意味であるので、「待ち行列」の先頭の機械カードを取り出し、修理者１のところにおく。そして、記録用紙の修理者１の欄に、修理する機械番号を１と書いておく。
⑤ 修理完了
もしも修理者１の位置に、機械カードが置いてあれば、それは修理者が修理中であること、つまり稼働中があることを意味している。「修理完了機構」のカー
ドをよくきってその中から１枚を引く。もし「修理完了」カードを引いたら、その時点で修理完了であるから、「修理者」の位置にあった機械カードをもとの
「機械」の位置にもどす。記録用紙の機械の欄に○（修理完了のマーク）をつける。
⑥ 修理完了繰返し
修理者２について④⑤を繰り返す。
⑦ 以上の事が、１時間の中で同時に起った事です。時計を１時間進め、②から繰り返す。
以上により、数十時間分を繰り返し、記録用紙を見れば、表８．１．１のようになるでしょう。各機械が何時間故障であったか、修理開始まで何時間待たされるか、修理者はどのくらい稼働していたか、この表から計算することができるでしょう。
この方法は、３．１節で述べたベルヌーイ試行の実験の応用で、時計を一時停めて、すべての機械、修理者に対し、ベルヌーイ試行を行った後、時計を進めると考えて下さい。これが同時進行のシミュレーション方法です。
| 図８．１．１ 机上実験 |
 |
| 表８．１．１ 記録用紙 |
 |
コンピュータで機械修理問題を解く
８．１節の机上実験で、同時進行事象のシミュレーションの概念はつかめたと思いますが、これを何百時間分も続けないと、信頼のおける機械の稼働率、修理
者の稼働率などは求められません。また、時刻の進め方も、１時間を単位にしていますが、３．３節で述べたように十分に微少な時間として、その中に２つの事
象の起こらないようにしなければなりません。このような目的のためには、コンピュータの手をかりるのが良いでしょう。そこで、同じ例題をコンピュータでシ
ミュレーションすることにしましょう。
プログラムs0820.cでは、次の変数を用います。机上実験で、机の上に置いていた「機械」カードは、 sm[j] に相当し、この値が０であれば、稼働中、１であれば故障中であることを示します。修理者の状態を
sn[k] で表わします。この値が、０であれば未稼働（ひまな状態）、１～５であれば機械１から機械５を修理していることを示します。この様子を図８．２．１に示します。このような図を状態遷移図と言います。また、
q は待ち行列の長さを、 sq[q] は待ち行列の q 番目にならんでいる機械の番号を表わすこととします。
| 図８．２．１ 機械修理問題の状態遷移図 |
 |
プログラムのポイント
プログラムs0820.cが、機械修理の問題をシミュレーションするプログラムです。８．１節の机上実験と対比すると、行34から行42までが、②、③に対応します。行36が「故障発生機構」に相当します。行38,39で待ち行列につけます。行45から行62までが８．１節の④、⑤、⑥に相当します。行46から行55が④の修理開始です。
待ち行列 sq[q] は図８．２．２のような構造となっており、q の値が待ち行列の最後尾を示していますから、 q+1 に機械番号を入れることは、待ち行列に新たにならぶことを意味し、 sq[1] をとり出せば待ち行列の先頭がとり出せます。行51から行53において、行列が１歩ずつ前へ進むことを示しています。行57は「修理完了機構」を示しています。行65以降は記録用紙への記入に相当し，行81は時刻を DT だけ進めます。
| 図８．２．２ 待ち行列の構造 |
 |
実行
実行すると、m=5, n=2, a=10, h=4の場合、図８．２．３の結果が得られ、DT=0.5としていますから、０．５時間ずつの記録がとられ、表８．１．１と対比することができます。
シミュレーションの最後にp[f]として、機械がf台故障している時間の割合が示されます。
| 図８．２．３ 固定時間方式シミュレーション結果(部分） |
 |
| 途中省略 |
 |
| s0820.c 機械修理の問題 (固定時間方式) |
/* s0820.c
* 機械修理の問題 (固定時間方式)
* (C) H.Ishikawa 1994 2018
*/
#include <stdio.h>
#include "compati.h"
#define M 5 /* 機械の台数 */
#define N 2 /* 修理者の人数 */
#define A 10.0 /* 平均故障発生間隔 */
#define H 4.0 /* 平均修理時間 */
#define DT 0.5 /* 時刻増分 */
#define T_END 10000.0 /* シミュレーションを終る時刻 */
int main(void)
{
double tt = 0.0; /* シミュレーションクロック */
int i; /* 待ち行列のforのカウンタ */
int j; /* 機械のforのカウンタ */
int k; /* 修理者のforのカウンタ */
int f = 0; /* 故障している機械の台数 */
int q = 0; /* 待ち行列の長さ */
int sm[M + 1] = {0}; /* 機械jの状態 */
int sn[N + 1] = {0}; /* 修理者kの状態 */
int sq[M + 1] = {0}; /* 待ち行列 */
double p[M + 1] = {0.0}; /* 機械がf台故障している時間 */
srand(15);
/* 開始 */
while (tt <= T_END) {
/* 故障発生 */
for (j = 1; j <= M; j ++) {
if (sm[j] == 0) { /* 機械が稼働中なら */
if ((1.0 / A * DT) > RAND()) {
sm[j] = 1; /* ある確率で故障する */
q = q + 1; /* 待ち行列の後ろにつく */
sq[q] = j;
}
}
}
/* 故障修理 */
for (k = 1; k <= N; k ++) {
if (sn[k] == 0) { /* 修理者が未稼働なら */
if (q > 0) {
sn[k] = sq[1]; /* 待ち行列から取り出す */
q = q - 1;
if (q > 0) {
for (i = 1; i <= q; i ++) {
sq[i] = sq[i + 1]; /* 待ち行列を1進める */
}
}
}
} else {
if ((1.0 / H * DT) > RAND()) { /* 修理者は機械を修理し */
sm[sn[k]] = 0; /* 稼働状態にする */
sn[k] = 0;
}
}
}
/* 表示 */
printf("%7.2f ",tt);
for (j = 1; j <= M; j ++) printf(" %2d",sm[j]);
printf(" ");
for (k = 1; k <= N; k ++) printf(" %2d",sn[k]);
printf(" ");
for (i = 1; i <= q; i ++) printf(" %2d",sq[i]);
printf("\n");
/* 故障している機械の台数を数える */
f = 0;
for (j = 1; j <= M; j ++) {
if (sm[j] == 1) { f = f + 1; }
}
p[f] = p[f] + DT;
/* 時刻を進める */
tt =tt + DT;
}
/* 結果出力 */
for (j = 0; j <= M; j ++) {
printf("p[%1d] =%10f\n", j, (double)p[j] / T_END);
}
return 0;
}
|
別な方法で機械修理問題を解く
８．２節のプログラムは、簡便ですが、次の点に難があります。それは、時刻
tt の進め方についてです。tt は、キザミ DT を平均故障間隔 a あるいは平均修理時間 h に比し十分小さくする必要があります。３．２節で述べたとおり、dt が大きいということは、ある時間内に発生する事象は２項分布となることを意味し、ランダムなポアソン分布とは異なってしまいます。
それでは DT を小さくするとすれば、何の変化もない時（故障も発生せず、修理完了もない時）にも、計算をしなければならず、シミュレーション時間が長くなります。
これに比べ、ここで述べる変化点方式は、図８．３．１ｂのように変化が起こるまでの時間を計算して、その間の時刻をスキップする方式で、何らかの状態変 化があるときだけ、状態変化のプログラムが動くようにします。そのためには、現在からあと何秒（何日）後に次の状態が変化するかという、時間に相当する乱 数を用います。これには４章で述べた方法が用いられます。例えば、先の例題で、ある機械が修理完了となったとき、次に故障するまでの時間は指数分布となり ますから、プログラムs0420.cで作った指数分布乱数を用い、次に故障するまでの時間を計算します。
ランダム発生の指数分布のみならず、４章で述べた各種の分布乱数を用いれば、さまざまな時間間隔もシミュレーションできることも、この方式の特徴です。
| 図８．３．１ 固定時間方式と変化点方式 |
 |
プログラムのポイント
さてプログラムs0830.cにより、変化点方式のシミュレーションの説明をします。プログラムs0820.cに比べ、プログラムがやや長くなっています。sm[j], sn[k], q, sq[q] は８．２節の場合と同じです。これらの変数に加えて、tm[j], tn[k] という変数をあらたに導入します。tm[j] は、機械jが sm[j] の状態にあと何時間とどまっているか、その残り時間を表わします。例えば、sm[1]=0 で、tm[1]=3.68 であれば、「機械１が現在稼働中で、あと3.68時間後に故障が発生する」ことを意味します。また
sn[2]=3 で tn[2]=1.20 であれば、「修理者2は、現在機械３を保守しており、あと１．２０時間後に修理完了する」ことを意味します。この
tm[j] あるいは tn[k] において、機械の修理を完了した、あるいは保守者が修理を開始するたびに指数乱数を発生させ、次に起こる変化点までの時間をセットすることとします。
つぎに時刻ttの進め方に関しては、すべての機械の tm[j] 、あるいは、すべての保守者の tn[k] の中から最も小さい値（これが最も近く起こる事象です）を探し、その値が０になるまで進めます。
行36から41で機械 1 から M を「稼働中」に、修理者 1 から N までを「未稼働」に初期設定します。ここで修理者の待ち時間を無限大に設定します。 行44からがメイン・プログラムで、機械1からMに対し状態の変化をさせます。行48で tm[j] が０の機械、すなわち故障した機械を見つけ行49で sm[j]=1 に変化させ、待ち行列にならべます。ここで、 tm[j] の値は浮動小数点ですから、何度も演算を繰り返している間に、完全に０でないことが想定されますので、 tm[j] の絶対値が EPS という十分小さい値より小さい場合を０とみなしています。
行56から行73が修理者の動きを記述しています。行57で tn[k]=0 の修理を終えた修理者を探し、行58でその修理していた機械を稼働状態にします。つぎに k の修理者がたった今修理を完了したか、以前から未稼働の場合は待ち行列を調べます。待ち行列に機械があれば、その中にならんでいる機械を取り出して保守します。行57から行72をすべての修理者に対し実行します。
行85から行91で、 tm[j], tn[k] の中で、最も小さい値を探し、それをatとします。この値が、この系の中で、次におこる状態変化までの時間となります。行101以下で、すべての tm[j], tn[k] から at を減じ、またクロック tt を at だけ進めます。
行76以下が途中経過を記録する部分、行94以下が故障している機械の台数を数える部分です。
実行結果
実行すると、m=5, n=2, a=10, h=4 の場合、図８．３．２の結果が得られます。時刻は等間隔ではなく進行しており、１行出力するたびに機械か、修理者の状態が変化しているのがわかります。
| 図８．３．２ 変化点方式シミュレーション結果 |
 |
| 途中省略 |
 |
| s0830.c 機械修理の問題 (変化点方式) |
/* s0830.c
* 機械修理の問題 (変化点方式)
* (C) H.Ishikawa 1994 2018
*/
#include <stdio.h> /* printf() */
#include <math.h> /* log() fabs() */
#include "compati.h"
#define M 5 /* 機械の台数 */
#define N 2 /* 修理者の人数 */
#define A 10.0 /* 平均故障発生間隔 */
#define H 4.0 /* 平均修理時間 */
#define T_END 10000.0 /* シミュレーションを終る時刻 */
#define INF 1.0e100 /* 十分大きい値 */
#define EPS 1.0e-100 /* 十分小さい値 */
int main(void)
{
double tt = 0.0; /* シミュレーションクロック */
int i; /* 待ち行列のforのカウンタ */
int j; /* 機械のforのカウンタ */
int k; /* 修理者のforのカウンタ */
int f = 0; /* 故障している機械の台数 */
int q = 0; /* 待ち行列の長さ */
int sm[10] = {0}; /* 機械jの状態 */
double tm[10]; /* 機械に次の状態変化が起るまでの時間 */
int sn[10] = {0}; /* 修理者kの状態 */
double tn[10]; /* 修理者に次の状態変化が起るまでの時間 */
int sq[10] = {0}; /* 待ち行列 */
double p[10] = {0.0}; /* 機械がf台故障している時間 */
double at; /* 時刻増分 */
/* 初期設定 */
srand(15);
for (j = 1; j <= M; j ++) {
sm[j] = 0; tm[j] = -A * log(1 - RAND());
}
for (k = 1; k <= N; k ++) {
sn[k] = 0; tn[k] = INF;
}
/* 開始 */
while (tt <= T_END) {
/* 故障発生 */
for (j = 1; j <= M; j ++) {
if (fabs(tm[j]) < EPS) { /* 機械jが故障発生 */
sm[j] = 1; tm[j] = INF;
q = q + 1;
sq[q] = j;
}
}
/* 故障修理 */
for (k = 1; k <= N; k ++) {
if (fabs(tn[k]) < EPS) { /* 修理者kの修理が完了 */
sm[sn[k]] = 0; tm[sn[k]] = -A * log(1 - RAND());
}
if (fabs(tn[k]) <EPS || sn[k] == 0) { /* kが未稼動の場合 */
if (q > 0) { /* 待ち行列を取り出す */
sn[k] = sq[1]; tn[k] = -H * log(1 - RAND());
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;
}
}
}
/* 表示 */
printf("%6.2f ",tt);
for (j = 1; j <= M; j ++) printf(" %2d", sm[j]);
printf(" ");
for (k = 1; k <= N; k ++) printf(" %2d", sn[k]);
printf(" ");
for (i = 1; i <= q; i ++) printf(" %2d", sq[i]);
printf("\n");
/* 次に状態変化の起る時刻を求める */
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];
}
/* 故障している機械の台数を数える */
f = 0;
for (j = 1; j <= M; j ++) {
if (sm[j] == 1) { f ++; }
}
p[f] = p[f] + at;
/* 時間を進める */
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 ++) {
printf("p[%1d] =%10f\n", j, (double)p[j] / T_END);
}
return 0;
}
|
固定時間方式と変化点方式
同じ例題を２つの方式（固定時間方式と変化点方式）でシミュレーションを試みました。これらの比較をしたいと思います。比較項目は計算機時間と、理論値に対する誤差です。 T_END=10000 時間分のシミュレーションを固定時間方式についてはdt=1 時間、0.1 時間、0.01 時間として実行し、変化点方式と比較することとします。
理論値
理論値は故障発生間隔、修理時間とも指数分布の場合、次のように求めることができます。機械がj 台こわれている状態をEj と書くと、状態遷移図は、図８．４．１のように書くことができます。ここでλ は単位時間あたりの故障確率、μ は単位時間あたりの修理完了確率で、λ=1/a、μ=1/h
です。この図で左から右へ移行するということは、機械がつぎつぎ故障するということで、その遷移は故障していない機械の台数に比例しておこります。従っ
て、 dt という微少時間の間に (m-j)λdt という確率で遷移することになります。右から左への遷移は、故障していた機械が、修理され回復していくという意味で、その割合は、j
が n より大きい場合は、nμdt という確率で、j が n より小さい場合は、現在修理中の台数に比例した jμdt という確率で修理完了します。
| 図８．４．１ |
 |
ある時刻 t で Ej という状態にある確率を pj(t) と書けば、この図から、微小時間経過後の pj(t+dt) を次のように、あらわすことができます。
 |
・・・ |
式８．４．１ |
ここで、定常状態の解を求めるため、
 |
・・・・・ |
式８．４．２式 |
とし、これを改めてpjと書くと、これらの式８．４．１は
 |
・・・ |
式８．４．３ |
を得ます。これらから、p0 からp1 ,p1 からp2 ,p3 からp4 というように求めていくと、
 |
・・・ |
式８．４．４ |
ただし、(p1 + p2 + p3 +・・・+ pm = 1)
という漸化式を得ることができます。
この式からプログラムs0840.cを書くことができます。このプログラムでは pj を p[j] 、m を M 、n を N 、1/λ を A 、1/μ を H と表わしています。
両方式の比較
前節で使った数値例でs0840.cにより、理論値を計算し、固定時間方式と変化点方式の比較をおこないます。図８．４．２が得ら れます。比較表を表８．４．１に示します。変化点方式は実行速度が1桁以上速いうえ、精度が良いことがわかります。また固定時間方式は、あとで述べる、在庫管理問
題のように１日単位に起こる変化をシミュレーションする場合などには、プログラムが単純となりますから、良く用いられます。問題によって使い分けるのがよ
ろしい。
| 図８．４．２ s0840.cの結果 |
 |
| 表８．４．１ 固定時間方式と変化点方式の比較 |
 |
| *.java は*.cと読み替えてください |
| s0840.c 機械修理の問題 (理論値) |
/* s0840.c
* 機械修理の問題 (理論値)
* (C) H.Ishikawa 1994 2018
*/
#include <stdio.h>
#define M 5 /* 機械の台数 */
#define N 2 /* 修理者の人数 */
#define A 10.0 /* 平均故障発生間隔 */
#define H 4.0 /* 平均修理時間 */
int main(void)
{
int j;
double p[10] ;
double p0 = 0.0;
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 ++) {
printf(" p[%1d] =%10f \n", j, p[j]);
}
return 0;
}
|
シミュレーションをおこなう場合，信頼に足る結果を得るのにどのくらい計算機を回すのが必要でしょうか。パソコンを用いれば，最終結果を アウトプットするのみならず，途中の経過をディスプレイすることができますから，その形をみて，収束状況を，判定することができます。
ディスプレイの画面には，横軸に試行回数あるいはシミュレーションクロック，縦軸に結果をとると，図８．５．１のようなグラフが表示されるでしょう。す
なわち，始めはばらつきが大きく，しだいに，収束されて行きます。３．４節の大数の法則と，３．５節の中心極限定理を，思いだして下さい。縦軸はある試行
の結果，成功の回数と試行の回数の比ですから，大数の法則により，試行回数
n
を多くすると一定値に近づくわけです。しからば，この場合のバラツキはどうでしょうか。中心極限定理は， x の分布にかかわらず， x を n 個合計
し， n で割ると，標準偏差が， １／√n の正規分布となるといっています。 n を大きくすると標準偏差は小さく中心に近ずくわけです。しかし，そ
の割合は n を，１００倍しても，標準偏差は，１／１０にしかなりません。つまり，１桁精度を上げようとすれば，シミュレーション時間を１００倍かける
必要があります。
| 図８．５．１ シミュレーション時間と結果の収束 |
 |
１． 機械の台数＝５台、故障の平均発生間隔＝１０時間、平均修理時間＝４時間、１時間あたりの機械１台の収益＝８万円、修理者１人１時間あたりの費用＝２万円としたとき、最大利益が得られる理論的な修理者の人数を求めなさい。
（解答） これによれば、２人が最良
| a0810.c ８章１．の解答 機械修理の問題 (理論値により最良の人数を求める ) |
/* a0810.c ８章１．の解答
* 機械修理の問題 (理論値により最良の人数を求める )
* (C) H.Ishikawa 1994 2018
*/
#include <stdio.h>
#define M 5 /* 機械の台数 */
#define A 10.0 /* 平均故障発生間隔 */
#define H 4.0 /* 平均修理時間 */
#define PROF 8.0 /* 機械１台当りの利益 */
#define COST 2.0 /* 修理者１組当りの費用 */
int main(void)
{
int j;
double p[10] ;
double p0 ;
double pr; /* 総利益 */
int n; /* 修理者の人数 */
for (n = 1; n <= M; n ++) {
pr = 0.0;
p0 = 0.0;
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 ++) {
pr = pr + p[j] * (double)(M - j) * PROF;
}
printf(" %3d %6.3f\n", n, pr - n * COST);
}
return 0;
}
|
２． プログラムs0830.cで、修理時間を位相k=3のアーラン分布とした場合について、シミュレーションしなさい。
（解答）
| a0820.c ８章２．の解答 機械修理の問題 (変化点方式) |
/* a0820.c ８章２．の解答
* 機械修理の問題 (変化点方式)
* (C) H.Ishikawa 1994 2018
*/
#include <stdio.h> /* printf() */
#include <math.h> /* log() fabs() */
#include "compati.h"
#define M 5 /* 機械の台数 */
#define N 2 /* 修理者の人数 */
#define A 10.0 /* 平均故障発生間隔 */
#define H 4.0 /* 平均修理時間 */
#define T_END 10000.0 /* シミュレーションを終る時刻 */
#define INF 1.0e100 /* 十分大きい値 */
#define EPS 1.0e-100 /* 十分小さい値 */
double erlang(double lambda, int k);
int main(void)
{
double tt = 0.0; /* シミュレーションクロック */
int i; /* 待ち行列のforのカウンタ */
int j; /* 機械のforのカウンタ */
int k; /* 修理者のforのカウンタ */
int f = 0; /* 故障している機械の台数 */
int q = 0; /* 待ち行列の長さ */
int sm[10] = {0}; /* 機械jの状態 */
double tm[10]; /* 機械に次の状態変化が起るまでの時間 */
int sn[10] = {0}; /* 修理者kの状態 */
double tn[10]; /* 修理者に次の状態変化が起るまでの時間 */
int sq[10] = {0}; /* 待ち行列 */
double p[10] = {0.0}; /* 機械がf台故障している時間 */
double at; /* 時刻増分 */
/* 初期設定 */
srand(15);
for (j = 1; j <= M; j ++) {
sm[j] = 0; tm[j] = -A * log(1 - RAND());
}
for (k = 1; k <= N; k ++) {
sn[k] = 0; tn[k] = INF;
}
/* 開始 */
while (tt <= T_END) {
/* 故障発生 */
for (j = 1; j <= M; j ++) {
if (fabs(tm[j]) < EPS) { /* 機械jが故障発生 */
sm[j] = 1; tm[j] = INF;
q = q + 1;
sq[q] = j;
}
}
/* 故障修理 */
for (k = 1; k <= N; k ++) {
if (fabs(tn[k]) < EPS) { /* 修理者kの修理が完了 */
sm[sn[k]] = 0; tm[sn[k]] = -A * log(1 - RAND());
}
if (fabs(tn[k]) <EPS || sn[k] == 0) { /* kが未稼動の場合 */
if (q > 0) { /* 待ち行列を取り出す */
sn[k] = sq[1]; tn[k] = erlang(1 / H, 3); /* アーラン分布 */
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;
}
}
}
/* 表示 */
printf("%6.2f ",tt);
for (j = 1; j <= M; j ++) printf(" %2d", sm[j]);
printf(" ");
for (k = 1; k <= N; k ++) printf(" %2d", sn[k]);
printf(" ");
for (i = 1; i <= q; i ++) printf(" %2d", sq[i]);
printf("\n");
/* 次に状態変化の起る時刻を求める */
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];
}
/* 故障している機械の台数を数える */
f = 0;
for (j = 1; j <= M; j ++) {
if (sm[j] == 1) { f ++; }
}
p[f] = p[f] + at;
/* 時間を進める */
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 ++) {
printf("p[%1d] =%10f\n", j, (double)p[j] / T_END);
}
return 0;
}
double erlang(double lambda, int k)
{
double tp = 1.0;
double tau;
int n;
for (n = 1; n <= k; n ++) {
tp = tp * (1 - RAND());
}
tau = -1.0 / lambda / (double)k * log(tp);
return (tau);
}
|
３． プログラムs0830.cを改造して、故障している機械の割合 pf[j]/tt を、時々刻々表示し、シミュレーションの収束状況を、観測できるようにしなさい。また機械の稼働率、修理者の稼働率、機械の待ち時間の分布、機械の待ち合わせ平均時間、その分散を求めるようにしなさい。
（解答）
| a0830.c ８章３．解答 機械修理の問題 (変化点方式 統計出力) |
/* a0830.c ８章３．解答
* 機械修理の問題 (変化点方式 統計出力)
* (C) H.Ishikawa 1994 2018
*/
#include <stdio.h> /* printf() */
#include <math.h> /* log() fabs() */
#include "compati.h"
#include "window.h"
#define M 5 /* 機械の台数 */
#define N 2 /* 修理者の人数 */
#define A 10.0 /* 平均故障発生間隔 */
#define H 4.0 /* 平均修理時間 */
#define T_END 10000.0 /* シミュレーションを終る時刻 */
#define INF 1.0e100 /* 十分大きい値 */
#define EPS 1.0e-100 /* 十分小さい値 */
int main(void)
{
double tt = 0.0; /* シミュレーションクロック */
int i; /* 待ち行列のforのカウンタ */
int j; /* 機械のforのカウンタ */
int k; /* 修理者のforのカウンタ */
int f = 0; /* 故障している機械の台数 */
int q = 0; /* 待ち行列の長さ */
int sm[10] = {0}; /* 機械jの状態 */
double tm[10]; /* 機械に次の状態変化が起るまでの時間 */
int sn[10] = {0}; /* 修理者kの状態 */
double tn[10]; /* 修理者に次の状態変化が起るまでの時間 */
int sq[10] = {0}; /* 待ち行列 */
double at; /* 時刻増分 */
double am[10] = {0}; /* 機械の稼働時間 */
double an[10] = {0}; /* 修理者の稼働時間 */
double tq[10] = {0}; /* 待行列にならんだ時刻 */
double pf[10] = {0}; /* 機械の故障している時間 */
int fr[50] = {0}; /* uの頻度 */
int ff = 0; /* ある時刻における故障機械の台数 */
int nn = 0; /* 修理した機械の台数 */
int u; /* twの整数部分 */
double tw; /* 待ち時間 */
double x = 0; /* twの和 */
double x2 = 0; /* twの自乗和 */
/* グラフィクの準備 */
opengraph();
SetWindow(0, 0.0,0.0,T_END,0.5, SPACE,HIGHT-SPACE,WIDTH-SPACE,SPACE);
Axis(0, "tt", T_END/10, "p[]", 0.1);
MoveTo(0, 0.0, 0.0);
/* 初期設定 */
srand(15);
for (j = 1; j <= M; j ++) {
sm[j] = 0; tm[j] = -A * log(1 - RAND());
}
for (k = 1; k <= N; k ++) {
sn[k] = 0; tn[k] = INF;
}
/* 開始 */
while (tt <= T_END) {
/* 故障発生 */
for (j = 1; j <= M; j ++) {
if (fabs(tm[j]) < EPS) { /* 機械jが故障発生 */
sm[j] = 1; tm[j] = INF;
q = q + 1;
sq[q] = j;
tq[q] = tt;
}
}
/* 故障修理 */
for (k = 1; k <= N; k ++) {
if (fabs(tn[k]) < EPS) { /* 修理者kの修理が完了 */
sm[sn[k]] = 0; tm[sn[k]] = -A * log(1 - RAND());
}
if (fabs(tn[k]) <EPS || sn[k] == 0) { /* kが未稼動の場合 */
if (q > 0) { /* 待ち行列を取り出す */
sn[k] = sq[1]; tn[k] = -H * log(1 - RAND());
tw = tt - tq[1]; u = (int)tw;
x = x + tw;
x2 = x2 + tw * tw;
fr[u] = fr[u] + 1;
nn = nn + 1;
q = q - 1;
if (q > 0) {
for (i = 1; i <= q; i ++) {
sq[i] = sq[i + 1];
tq[i] = tq[i + 1];
}
}
} else { /* 待ち行列が空なら修理者kは未稼働状態になる */
sn[k] = 0; tn[k] = INF;
}
}
}
/* 次に状態変化の起る時刻を求める */
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];
}
for (j = 1; j <= M; j ++) {
tm[j] = tm[j] - at;
if (sm[j] == 1) { ff = ff + 1;}
if (sm[j] == 0) { am[j] = am[j] + at;}
}
for (k = 1; k <= N; k++) {
tn[k] = tn[k] - at;
if (sn[k] > 0) { an[k] = an[k] + at;}
}
pf[ff] = pf[ff] + at; ff = 0;
tt = tt+ at;
/* 表示 */
/* printf("%6.2f ",tt);
for (j = 1; j <= M; j ++) printf(" %2d", sm[j]);
printf(" ");
for (k = 1; k <= N; k ++) printf(" %2d", sn[k]);
printf(" ");
for (i = 1; i <= q; i ++) printf(" %2d", sq[i]);
printf("\n");
*/
/* グラフ表示 */
for (j = 0; j <= M; j ++) {
PutPixel(0, tt, ((double)pf[j] / tt), (j + 1));
}
}
/* 出力 */
printf("故障している時間の割合 \n");
for (j = 0; j <= M; j ++) {
printf(" p[%2d] =%10f \n", j, (double)pf[j] / T_END);
}
printf("機械稼働率 \n");
for (j = 1; j <= M; j ++) {
printf(" am[%2d] =%10f \n", j, (double)am[j] / T_END);
}
printf("修理者稼働率 \n");
for (k = 1; k <= N; k ++) {
printf(" an[%2d] =%10f \n", k, (double)an[k] / T_END);
}
printf("待ち時間の分布\n");
for (u = 0; u <= 9; u ++) {
printf(" fr[%2d] =%10f \n", u, (double)fr[u] / nn);
}
printf("\n");
printf("平均待ち時間 =%10f\n", x / nn);
printf("待ち時間の分散 =%10f\n", (x2 - x * x / nn) / nn);
getch();
closegraph();
return 0;
}
|
４． １本の動力線を m 人の溶接工が共用しています。溶接工はランダムに断続的に１Ａの電流を使用します。その時間の平均値は h 分で間隔の平均は a 分です。電流の波形をシミュレーションしなさい。
（解答）
| a0840.c ８章４．解答 溶接の電流 |
/* a0840.c ８章４．解答
* 溶接の電流
* (C) H.Ishikawa 1994 2018
*/
#include <math.h>
#include "compati.h"
#include "window.h"
#define M 20
#define H 10
#define A 20
#define T_END 1000 /* 終りの時刻 */
#define INF 1.0e100 /* 十分大きい値 */
#define EPS 1.0e-100 /* 十分小さい値 */
int main(void)
{
double tt = 0.0;
int j;
int s[M + 1];
double t[M + 1];
int mm = 0; /* 同時使用数（電流値） */
double at; /* 次に状態の変化が起るまでの時間 */
double cp = 0.0; /* カレントポインタ */
/* グラフィクの準備 */
opengraph();
SetWindow(0, 0.0,0.0,T_END,M, SPACE,HIGHT-SPACE,WIDTH-SPACE,SPACE);
Axis(0, "tt", T_END/10, "mm", 1);
MoveTo(0, 0.0, 0.0);
/* 初期設定 */
srand(3);
for (j = 1; j <= M; j ++) {
s[j] = 0; t[j] = -A * log(1.0 - RAND());
}
/* メイン */
while (tt <= T_END) {
for (j = 1; j <= M; j ++) {
if (fabs(t[j]) < EPS) {
if (s[j] == 0) {
s[j] = 1; t[j] = -H * log(1.0 - RAND());
mm = mm + 1;
} else {
s[j] = 0; t[j] = -A * log(1.0 - RAND());
mm = mm - 1;
}
}
}
LineTo(0, tt, cp);
LineTo(0, tt, (double)mm);
cp = (double)mm;
/* 時間を進める */
at = INF;
for (j = 1; j <= M; j ++) {
if (t[j] < at) at = t[j];
}
for (j = 1; j <= M; j ++) {
t[j] = t[j] - at;
}
tt = tt+ at;
}
getch();
closegraph();
return 0;
}
|
|
 |
