 |
本書の最後は、「何でもシミュレーションしてみよう」と題し、ギャンブルから、最近注目されている、ジェネティックアルゴリズムまで、乱数を用いた、さまざまな、確率実験を、試みます。
クラップスという名のゲーム文献（８）（１３）
カジノで行われているダイス（サイコロ）を使った、クラップスというゲームをとりあげます。このクラップスに限らず、カジノで行われるゲームは、すべて
胴元が有利となるようなルールが定められています。したがって長時間平均すれば、必ず胴元が勝つ仕組みとなっているわけですが、必勝法はないものでしょう
か。
まず、クラップスのルールを説明します。クラップスは２個のサイコロを用いておこないます。２個のサイコロの目の合計は２から１２までの値をとりま
す。（３．５節で示したとおり、７の出る確率が最も大きい。）プレーヤは２個のサイコロを振って、その和が７または、１１の時（これをナチュラルといいま
す）は勝ちで,２、３あるいは１２の時（これをクラップスといいます）は負けです。もしそれ以外のとき、すなわち、４、５、６、８、９または、１０がでた
ときは、その値をポイントといい、ひきつづき同じポイントまたは７が出るまで、サイコロを振ります。そして同じポイントがでれば、プレーヤの勝ち、それ以
前に７が出れば、負けとなります。もしプレーヤが勝てば、自分の賭け金の２倍の払い戻しを受けます。
このゲームは、たいへんよくできて、ごく僅かな確率の差で、胴元が勝つようにできています。
必勝法
このゲームで必ずプレーヤが勝つ方法はあるでしょうか。もしそれがあれば、ぜひカジノに行って見たいものです。倍増法という必勝法が知られています。倍
増法というのは、まず１ドルの賭け金からスタートし、もし勝てば、次も１ドルを賭けますが、負けた場合は２ドルを、さらに負ければ前回の賭け金の２倍の額
を次々に賭けていく方法です。この方法では、負けた金額をとりもどせるので、必勝法となりそうです。はたして、これがほんとうに必勝法となるのでしょう
か。クラップスにおいて、倍増法をシミュレーションしてみましょう。
シミュレーション・プログラム
ではクラップスのシミュレーションプログラムs1010.cの説明を行ないます。まず、サイコロの関数を作ります。
int dice(void);
により、整数１～６の一様分布の乱数を発生させます。次に、(2)のとおりクラップスのルールどおりに記述して、関数
int craps(void);
を作ります。１回目のサイコロの目の和を sum1 , ２回目以後のサイコロの目の和を sum2 としています。この関数は、勝ちの場合１、負けの場合０の値を返します。
次にメインプログラムですが、１回の賭け金を k ドル、現在の持ち金を q ドルとします。行30から行35が倍増法による賭けで、もし勝った場合、次の賭け金を１とし、負けた場合前回の２倍を賭けることとします。プログラムでは、賭けを何回も繰り返し、横軸賭けの回数 t 、縦軸持ち金 q をグラフに表示するようにしています。
シミュレーションの結果
プログラムを実行すると、図１０．１．１のように時間とともに、着実に持ち金は増加していきますが、時々負けが続くと賭け金が急増し、赤字に転落すると
きがあります。もしその時、もち金がなくなれば破産ですが、なお余裕があって倍々と賭けることができれば、必ずプラスになります。しかしながら、通常は胴
元のほうが資金は豊富で、負けが続いても持ちこたえられるのに対し、客のほうが先に底をついてしまいます。もし十分な資金を用意することができれば、ぜひ
カジノへ行って、ひと儲けして下さい。
| 図１０．１．１ 倍増法による必勝法 |
 |
| s1010.c ギャンブル必勝法 |
/* s1010.c
* ギャンブル必勝法
* (C) H.Ishikawa 1994 2018
*/
#include <time.h> /* time() */
#include "compati.h"
#include "window.h"
#define T_END 500
long t; /* クロック */
long q = 100; /* 持ち金 */
long k = 1; /* １回の賭金 */
int dice(void);
int craps(void);
int main(void)
{
opengraph();
SetWindow(0, 0.0,0.0,(double)T_END,500.0,
SPACE,HIGHT-SPACE,WIDTH-SPACE,SPACE);
Axis(0, "t", T_END/10, "q", 100);
MoveTo(0, 0.0, q);
srand((unsigned)time(NULL)); /* シミュレーションの都度乱数系列を変える */
for (t = 1; t <= T_END; t ++) {
q = q - k;
if (craps() == 1) { /* 勝った場合 */
q = q + 2 * k;
k = 1; /* １ドル賭ける */
} else { /* 負けた場合 */
k = k * 2; /* ２倍賭ける */
}
LineTo(0, (double)(t), (double)(q));
}
getch();
closegraph();
return 0;
}
int dice(void)
{
return((int)(RAND() * 6.0 + 1));
}
int craps(void)
{
int r; /* 1:勝ち 0:負け */
int sum1; /* 和１ */
int sum2; /* 和２ */
sum1 = dice() + dice();
if (sum1 == 7 || sum1 == 11) { /* ナチュラル */
r = 1;
} else if(sum1 == 2 || sum1 == 3 || sum1 == 12) { /* クラップス */
r = 0;
} else {
sum2 = 0;
r = 0;
while (sum2 != 7) { /* ７が出るまで */
sum2 = dice() + dice();
if (sum2 == sum1) { /* ポイント */
r = 1;
break;
}
}
}
return(r);
}
|
在庫管理問題とは
オペレーションズリサーチの教科書に関らず出てくる問題です。企業や産業システムでは、原料や製品が供給から需要ひ向かって流れていく過程において、それらを貯蔵する機能が必要です。この貯蔵容量と、それを運用管理する最適の方法を決定するのが、在庫管理問題です。
例えば、ある商店が商品を問屋から仕入れて、客に販売する場合を考えてみましょう。商品を一度に多量仕入れると、在庫をたくさんかかえ、そのための費用
がかさみます。逆に少量の在庫では、大きな需要があった場合に品切れの危険が生じ、得べかりし利益をそこない、ひいてはその店の信用を失います。しから
ば、① 在庫に要する費用、② 問屋への発注のための費用、③ 品切れが発生した場合にこうむる損失の合計を最小にして、利益を最大にするには、いつどれ
だけ問屋に発注すべきでしょうか。
販売店の在庫 q は、図１０．２．１のように毎日の売上げ k 個ずつ減少していきます。基準発注点 q1 をわったら、 a 個の発注をしますが、すぐには入荷せず、 lt 日の遅れがあります。図１０．２．２のように、在庫量は変化するでしょう。売上げ個数と入荷おくれが確率的である場合を考え、 基準発注点q1 の最適値を求めるためのシミュレーション・プログラムを作ってみましょう。
| 図１０．２．１ 在庫モデル |
 |
| 図１０．２．２ 在庫 q の変化と lt |
 |
プログラムのポイント
売上げ個数は、毎日平均 AL のポアソン分布に従うとします。また入荷遅れ日数 lt は平均 ET 、標準偏差 ST の正規分布にしたがうものとします。これらの乱数の発生には４章で作った関数を用います。シミュレーション・プログラムの方式は１日をクロックとして進める固定時間方式を用います。プログラムs1020.cの行60から行78が，１日分を疑似します。ここで重要な変数は、入荷遅れ日数 lt で次の意味を持ちます。
lt > 0 ：すでに発注ずみであり、あと lt 日で納品される。
lt = 0 ：本日納品される。
lt < 0 ：発注していない。
行61で１日平均 k 個ポアソン分布の売上げがあったとします。行64で発注点 q1 を割り、まだ発注していなければ発注します。行68で q が負になれば、在庫なしで販売しそこなった数 ls を数えます。行73では、本日納品されたので、在庫に A を加えます。行75で在庫費用の累計を計算します。 T_END 日分のシミュレーションを終わったあと、発注費用の累計 sc2 、品切れ損失の累計 sc3 、売上個数の累計 yk を計算し、利益 rr を求めます。そしてグラフをプロットし、発注点 q1 をかえ、シミュレーションを続けます。
アウトプット
１００００日分のシミュレーションを q1 を０から１００まで変化させ、繰り返し行います。 発注点 q1に対する総売上げ rrのグラフ図１０．２．３を見ると、q1=45 の時に、利益が最大となることがわかります。
| 図１０．２．３ 在庫管理の問題 シミュレーション結果 |
 |
| s1020.c 在庫管理の問題 |
/* s1020.c
* 在庫管理の問題
* (C) H.Ishikawa 1994 2018
*/
#include <stdio.h>
#include <math.h>
#include "compati.h"
#include "window.h"
#define A 100 /* １回の入荷量 */
#define AL 5 /* １日の平均売上個数 */
#define ET 6 /* 納品するまでの平均日数 */
#define ST 2 /* 納品するまでの日数の標準偏差 */
#define SP 7000.0 /* 売値 */
#define PC 5000.0 /* 仕入れ単価 */
#define C1 20.0 /* １日１個当たりの在庫費用 */
#define C2 10000.0 /* １回の発注にともなう費用 */
#define C3 500.0 /* 品切れ１個あたりの損失 */
#define T_END 10000
#define Q1_MIN 0
#define Q1_MAX 90
long poisson(double lambda);
double normal(double ex, double sd);
int main(void)
{
long tt; /* シミュレーション日数 */
long q; /* 現在在庫量 */
long q1; /* 発注点 */
long k; /* １日の売上個数 */
long zk; /* kの累計 */
long yk; /* 売上個数の累計 */
long ls; /* 在庫割れのため売り損なった個数 */
long nn; /* 発注回数 */
long lt; /* 入荷遅れ日数 */
double sc1; /* 在庫費用の累計 */
double sc2; /* 発注費用の累計 */
double sc3; /* 品切れ損失の累計 */
double rr; /* 利益 */
/* グラフィック準備 */
opengraph();
SetWindow(0, 0.0,0.0,Q1_MAX,1.0e08,
SPACE,HIGHT-SPACE,WIDTH-SPACE,SPACE);
Axis(0, "q1", 10, "rr", 1e07);
setlinestyle(0, 0, 2);
//printf(" q1 yk ls nn sc1 sc2 sc3 rr\n");
/* 開始 */
for (q1 = Q1_MIN; q1 <= Q1_MAX; q1 = q1 + 5) {
/* 初期設定 */
srand(113);
tt = 0; q = A; zk = 0; ls = 0; nn = 0; lt = -1;
sc1 = 0.0; sc2 = 0.0; sc3 = 0.0;
/* １日分のシミュレーション */
while (tt <= T_END) {
k = poisson(AL); /* k個売れた */
q = q - k;
zk = zk + k;
if (q < q1 && lt < 0) { /* 発注点を割ったので発注する */
lt = (long)normal(ET, ST);
nn = nn + 1;
}
if (q < 0) { /* 在庫割れのため売り損ない */
ls = ls - q;
q = 0;
}
if (lt == 0) { /* 入荷 */
q = q + A;
}
sc1 = q * C1 + sc1;
lt = lt - 1;
tt = tt + 1;
}
sc2 = nn * C2;
sc3 = ls * C3;
yk = zk - ls;
rr = (SP - PC) * yk - sc1 - sc2 - sc3;
//printf("%4ld %6ld %6ld %6ld %8.0f %8.0f %8.0f %8.0f\n", q1, yk, ls, nn, sc1, sc2, sc3, rr);
if (q1 == Q1_MIN) {
MoveTo(0, q1, rr);
} else {
LineTo(0, q1, rr);
}
}
getch();
return 0;
}
long poisson(double lambda)
{
double xp = RAND();
long k = 0;
while (xp >= exp(-lambda)) {
xp = xp * RAND();
k = k + 1;
}
return (k);
}
double normal(double ex, double sd)
{
double xw = 0.0;
double x;
long n;
for (n = 1;n <= 12; n ++) { /* 12個の一様乱数の合計 */
xw = xw + RAND();
}
x = sd * (xw - 6.0) + ex;
return (x);
}
|
鉄は強磁性体
本節では、シミュレーション技術を用いて、１つ１つの原子の動きを計算し、物質全体としてどういう性質をもっているかを、鉄のような強磁性体を例に解析してみましょう。
原子の模型は、図１０．４．１のように、原子核のまわりを電子がまわっているように書けます。電子はマイナスの電気をもっており、電子が回転運動すると
磁場ができます。その軸は原子の中心を通り、電子の軌道面に対して垂直です。この磁場のことを電子の軌道磁気モーメントと言います。また電子は自転しなが
ら軸のまわりを回っており、これによっても、磁場が生じています。これをスピン磁気モーメントと言います。これら二種類の磁気モーメントが加え合わされ、
原子１つの磁気モーメントとなります。
| 図１０．４．１ 原子の模型 |
 |
鉄の原子の場合は、スピン磁気モーメントの釣合がひどくずれており、原子１つ１つが小さな磁石となっており、その強さはすべての元素の原子なかで一番強
いのです。鉄の原子は、結晶の格子点に固定されて動くことができませんが、その向きは自由に変り、いろいろの方向を向きます。鉄の棒を磁化して磁石にする
というのは、できるだけたくさんの原子の磁場の方向を並行にそろえることにほかなりません。
これはごく大ざっぱにいうと、鉄原子の電子スピンによる磁気モーメントは、近接する原子間では、互いに逆向きよりも平行の方がエネルギーが低くなるような相互作用が動いているためです。（図１０．４．２ａ）
磁石が磁石でなくなる
ところが温度を高くすると、原子はしだいに無秩序になり、上向き下向きのスピンの原子があらわれます。そしてある温度（ピエル・キュリーが発見した。
キュリー温度という）を越すと、向きがばらばらになり、鉄全体の磁気モーメントはほぼ零になり、鉄は磁石でなくなってしまうのです（図１０．４．２ｂ）。
このことを、強磁性体からの相転移といいます。シミュレーションにより、この相転移の実験をしてみましょう。
| 図１０．４．２ 絶対零度の場合とある温度の場合のスピン磁気モーメント |
 |
イジング・モデル
今原子が図１０．４．３のように２次元格子状に並んでいるとします。原子
(i,j) の状態は、スピンが上向きか、下向きかという２つの状態だけをとるものとし、それを
s(i, j)=1 , s(i, j)= -1 とします。この値をスピン変数といいます。
ある原子 (i0, j0) が、まわりの原子から受けるエネルギーは、まわりの原子のスピンの方向と、自分のスピンの方向によって決まり、その値は、
 |
・・・・・ |
式１０．４．１ |
となります。 J は物質によって異なり、 (i0, j0) から離れるにしたがって小さな値となります。特に原子 (i0, j0) から見てすぐ隣りの原子からのみ相互作用を受けるとすれば、式は、
 |
|
|
| ・・・式１０．４．２ |
|
|
となります。 e の値は、原子 (i0, j0) から見て、まわりの原子のスピンの方向が全部逆であれば、一番大きく、また原子
(i0, j0) とまわりの原子のスピンが全部同じであれば、一番小さくなります。
ある温度で、原子 (i0, j0) が向きを反転する確率は、統計力学の教えるところにより、単位時間あたり、
 |
・・・・・ |
式１０．４．３ |
となります。ここで、 k はボルツマン定数、 T は絶対温度、 e は式の値、 A は定数です。
式１０．４．２と、式１０．４．３から、原子 (i0, j0) が反転する確率は、 T が高くなると高くなり、逆向きの原子にとりかこまれているとき最も高く、まわりの原子と同じ向きのとき最も低くなります。
このようなモデルを提案者の名前をとって、イジング・モデルといいます。
| 図１０．４．３ |
 |
プログラムのポイント
以上のモデルを、シミュレーションにより実現します（プログラムs1040.c）。温度を高温にすると、あちらこちらで、同時に反転が起こりますが、シミュレーションでは、時計を停止させ、ランダムに原子を選び出し、処理します。
それにはまず、モデルの磁石の原子に相当する配列を用意します。行14で１辺の大きさが、 MAX+2 （値は 0 から MAX+1 ）の配列を作り、その中の 1 から MAX までを使います。 0 と MAX+1 は式１０．４．２の計算をするとき i 、 j が 1 や MAX の値をとったとき、その隣りの値が配列をはみださないようにするためです。そして初期状態は、絶対零度としすべて磁気モーメントは同方向、すなわち
s(i,j)=1 とします。
次に行47，行48で２つ１組の １ から ＭＡＸ までの整数値をとる一様分布の乱数を発生させます。この１組の乱数により、原子を１つ選び、式１０．４．２の相互作用のエネルギーを計算し（行49）、式１０．４．３の反転の確率を計算します（行50）。ここで
w0 = kT/J とおいて正規化温度としています。また A = exp(-4/w0) とします。別に発生した乱数とこの確率を比べ（行51）、乱数の方が小さければ、原子を反転させます（行52）。この処理を ＭＡＸ×ＭＡＸ 回繰り返すと、各原子について、平均１回処理を行ったことになるので、クロックを１進めます。
アウトプット
s(i, j)=1 である原子の数を数えて、磁化率を計算します（行62）。磁化率は全部の原子が s(i, j)=1 であれば１, １と－１が半々であれば０の値を取ります。また、ディスプレイに、刻々と原子の反転の様子をプロットさせます。それには反転のたびに
s(i, j) が１の時は青、－１の時は赤で表示します。これによりミクロな原子の動きをパターンとして視覚に訴えることができます。
プログラムを実行してみましょう。図１０．４．４のように絶対零度では、一面に正のモーメント（図では青）であったものが、しだいに負のモーメントの部分（図では赤）ができ、時間の経過とともにそれが合体し（これを磁気ドメインといいます）、形も刻々変化していきます。
ウィンドウ１のグラフは、磁化率をあらわします。
ある温度（キュリー点）以下では、磁化率は０にならず磁石のままで、ある温度
Tc をこえると、相転移が起こり磁化率が０となり、磁石が磁石でなくなっていきます。この時の正規化温度の、理論値は、
 |
・・・・・ |
式１０．４．４ |
で与えられます。 プログラムでは正規化温度w0 を変えることができるようになっていますから、 wc を中心に色々変えて相転移が起こる様子を実験してみて下さい。
| 図１０．４．４ イジングモデル シミュレーション結果 |
 |
| しだいに負のモーメントの部分が増加（図では赤） |
 |
| 時間の経過とともに、うごめうごめきながらそれが合体していく |
| s1040.c イジングモデル |
/* s1040.c
* イジングモデル
* (C) H.Ishikawa 1994 2018
*/
#include <math.h> /* exp() */
#include "compati.h"
#include "window.h"
#define MAX 120 /* 配列の大きさ */
#define T_END 1000 /* シミュレーションを終わる時刻 */
int s[MAX+2][MAX+2]; /* 磁気モーメントの値
s[i][j]=1あるいはs[i][j]=ｰ1 */
int i,j; /* 着目している原子の番号 */
int e; /* 相互作用の値 */
int m; /* s[i][j]=1である原子の数 */
int n; /* forのカウンタ */
int tt = 0; /* シミュレーションの時刻 */
double p; /* 反転の確率 */
double g; /* 磁化率 */
double w; /* 温度の相対値 */
int main(void)
{
/* 準備 */
opengraph();
SetWindow(0, 0,0,MAX,MAX, WIDTH/2,MAX*2,WIDTH/2+MAX*2,0);
Axis(0, "", MAX, "", MAX);
SetWindow(1, 0,0.0,T_END,1.0, SPACE,HIGHT-SPACE,WIDTH-SPACE,MAX*2+SPACE);
Axis(1, "tt", T_END / 10.0, "g", 0.1);
MoveTo(1, 0, 1.0);
w = 2.8; /*この値を変更して実験*/
m = MAX * MAX;
for (j = 0; j <= MAX ; j ++) {
for (i = 0; i <= MAX ; i ++) {
s[i][j] = 1;
}
}
/* メイン */
while (tt <= T_END) {
for (n = 1; n <= MAX * MAX; n ++) {
i = (int)(MAX * RAND()) + 1; /* 1からMAXまでの乱数 */
j = (int)(MAX * RAND()) + 1;
e = -s[i][j] * (s[i][j - 1] + s[i][j + 1] + s[i - 1][j] + s[i + 1][j]);
p = exp((e - 4.0) / w); /* 反転確率 */
if ( p > RAND()) {
s[i][j] = -s[i][j]; /* 反転 */
m = m + s[i][j];
}
if (s[i][j] == 1) {
PutPixel(0, i, j, 1); /* 青 */
} else {
PutPixel(0, i, j, 12); /* 明るい赤 */
}
}
tt = tt + 1;
g = (double)(2 * m) / (double)(MAX * MAX) - 1.0;
setcolor(15);
LineTo(1, tt, g);
}
getch();
closegraph();
return 0;
}
|
ジェネティック・アルゴリズムとは
ジェネティック・アルゴリズム（ＧＡ）は遺伝的アルゴリズムともよばれ、生物の進化過程をシミュレーションすることによって、たくさんの組み合せのなか
から最適な解を見つけ出す手法のことをいいます。ダーウィンの自然淘汰理論というのがありますが、これを基礎に、生物の進化過程を工学的モデルにおきかえ
たものです。
進化過程の概要は、おおよそ次のようにあらわされるでしょう。まず、０と１の値をもつ遺伝子(Gene)からなる染色体（Chromosome)を考え
ます。この染色体のあつまりを集団(Population)といいます。それに対し図１０．５．１ように（１）、（２）、（３）の操作をほどこします。
（１） 交叉
集団内の２つの染色体の組（ペア）において、交叉(crossover)という遺伝子の一部の入れ替えがおこなわれること。
（２） 突然変異
各染色体において、ある確率で、遺伝子の一部が０から１あるいは、１から０に返転すること。
（３） 選択
（１）、（２）をへた遺伝子をもった染色体のうち環境へのより高い適応性をもつものが、選択されること。
（１）、（２）、（３）の遺伝子的操作を１世代とよび、この操作をくりかえしていくと、その集団はより環境に順応する方向に進化します。
ジェネティック・アルゴリズムでは、さまざまな組み合わせを遺伝子からなる染色体で表現し、その組み合わせごとの目的関数を計算します。そして、その値
の高いものを選択し、それを残すよう生存競争を何世代にもわたり行わせれば、最終的には、目的関数の値が最大となるようにすることができます。
ジェネティック・アルゴリズムの原理は、驚くほど単純ですが、その探索能力は非常に高く、通常の方法では解を見つけにくい場合でも、比較的容易に正確に近づくことができるといわれています。
| 図１０．５．１ 組み替え（交叉と突然変異） |
 |
ナップザック問題
ジェネティック・アルゴリズムのメカニズムを「ナップザック問題」という典型的な組み合せ問題を例に、プログラムを作ってみましょう。
図１０．５．２のように、ナップザックがあります。ナップザックにはこべる重さには限界があります。たくさんの異なる重さの、異なる価値をもつ品物が与
えられたとき、重さがナップザックの限界値内で品物を幾つか選択し、その時の価値の合計が最大となるような組み合せを見つけるという問題です。すなわち、
この場合の目的関数は価値の合計で、ナップザックがこわれない範囲で、価値が最も高くなるよう、品物を選んで運ぼうとするものです。通常の方法では解を見つけにくい例です。
| 図１０．５．２ ナップザックの問題 |
 |
ナップザック問題のプログラム
ナップザック問題を、ジェネティック・アルゴリズムの手法で解く場合、まず、ナップザックに品物が入っている状態を染色体に対応させて、表現する必要があります。図１０．２．５の例では、品物１、３、４をナップザックに入れている状態ですから、これを染色体で
０１０１１
と表現します。
プログラムs1050.cについて説明していきます。LENGTHという長さをもつ染色体をSIZE個用意します。実際には組み替え操作のため，この２倍の配列を行19により用意します。また，行20により新世代の染色体を表わす配列も準備します。このほか１つ１つの品物の重量，価値をあらわす配列，染色体に対応した重量の合計，価値の合計をあらわす配列を用います。
行50から行56までで染色体を初期状態にします。PROB2の確率で１にセットし，また行59から行62において，各品物ごとの重量と価値も乱数でセットしておきます。
行65からがメインの始まりで，遺伝子の組み替え（交叉と突然変異），目的関数の計算，ソート，新世代へのおきかえ，表示という手続きを，G_ENDで示される世代にわたって繰り返し計算します。
遺伝子の組み替えは，図１０．５．３のように行います。k番目の染色体と，k+SIZE番目の染色体がペアであるとして，交叉の操作を行います。行69のように交叉位置を乱数で発生させ，k番目の後半とk+SIZE番目の後半を入れ換えます。次に行85以下のようにkの値がSIZE以上の染色体についてPROBという確率で，１と０を反転させ，突然変異を起こさせます。
次に，目的関数を計算します。目的関数は，染色体ごとの，合計の重量と合計の価値です。gene[k][l]に選ばれた品物の組み合せが示されていますから，行98あるいは行99のように，その配列に重量，あるいは価値の配列をかけ，合計の重量，価値を計算します。その場合，重量の合計が，ナップザックの重量の上限LIMITをこえた場合，行102により価値をゼロとし，その染色体が選択されないようにします。
次に，行107以下のように，価値の合計を大きい順にソートします。ソートのアルゴリズムはさまざまなものがありますが，ここでは，簡単な選択法というアルゴリズムを用います。total_value[k]について，まず，kを固定し，k+1以下のものと比較し，より大きいものがあれば，それと置き換えます。そして，kを次々増やして行くと，total_value[k]が大きい順に並びます。また，v_sort｢k]に，0,1,2,3,･･･の値を入れておくと，このソートが終わった状態でv_sort[k]に,total_value[k]の大きい方から順番の値が入ります。
ソートの結果を示すv_sort[k]により，行126のように，大きい順に，gene[k][l]を取り出し，それを新世代のnext_gene[k][l]の染色体を作り，行131によりgene[k][e]を入れ換えます。
| 図１０．５．３ 遺伝子の組み換え |
 |
結果の出力
２つのウインドウに横軸世代、縦軸それぞれ重量の合計、価値の合計をグラフ表示します。また、もっとも目的関数の大きい k = 0 の染色体の構造、またそのときの重量の合計、価値の合計を表示します。
このプログラムを走らせると、図１０．５．４のように、１００世代ほどで、最適に近い値が得られていますが、１０００世代くりかえし実行する必要があります。
| 図１０．５．４ ナップザック問題 シミュレーション結果 |
 |
| s1050.c ナップザック問題 |
/* s1050.c
* ナップザック問題
* (C) H.Ishikawa 1994 2018
*/
#include <stdio.h>
#include "compati.h"
#include "window.h"
#define LENGTH 40 /* 染色体の長さ */
#define SIZE 10 /* 集団の大きさ */
#define MAX_WEIGHT 10.0 /* 品物の重量の最大値 */
#define MAX_VALUE 10.0 /* 品物の価値の最大値 */
#define LIMIT 100.0 /* ナップザックに入る重量の限界値 */
#define G_END 2000 /* 繰り返し数 */
#define PROB 0.01 /* 突然変異の確率 */
#define PROB2 0.5 /* 初期化の確率 */
int gene[2*SIZE][LENGTH] = {0}; /* 染色体 */
int next_gene[2*SIZE][LENGTH] = {0}; /* 新世代の染色体 */
double weight[LENGTH]; /* 品物の重量 */
double value[LENGTH]; /* 品物の価値 */
double total_weight[2*SIZE]; /* 重量の合計 */
double total_value[2*SIZE]; /* 価値の合計 */
int v_sort[2*SIZE]; /* 順番 */
long g; /* 世代を表すforのカウンタ */
int k,l; /* gene[k][l]のためのforのカウンタ */
int l_rand; /* 交叉位置 */
int k_sort; /* ソート用のforのカウンタ */
int i_swap; /* スワップ用変数 */
double swap; /* スワップ用変数 */
int main(void)
{
/* グラフィックの準備 */
opengraph();
SetWindow(0, 0.0,0.0,G_END,120,
SPACE,HIGHT/2-SPACE/2,WIDTH-SPACE,SPACE);
Axis(0, "generation", G_END / 10, "weight", 20);
MoveTo(0, 0.0, 0.0);
SetWindow(1, 0.0,0.0,G_END,250,
SPACE,HIGHT-SPACE,WIDTH-SPACE,HIGHT/2+SPACE/2);
Axis(1, "generation", G_END / 10, "value", 50);
MoveTo(1, 0.0, 0.0);
setlinestyle(0, 0, 2);
srand(3);
/* 染色体のイニシャライズ */
for (k = 0 ; k < SIZE; k ++) {
for (l = 0; l < LENGTH; l ++) {
if (RAND() < PROB2){
gene[k][l] = 1;
}
}
}
/* 品物の重さと価値の設定 */
for (l = 0; l < LENGTH; l ++) {
weight[l] = MAX_WEIGHT * RAND();
value[l] = MAX_VALUE * RAND();
}
/* メイン始まり */
for (g = 1; g <= G_END; g ++) {
/* 交叉 */
for (k = 0; k < SIZE; k = k + 2) {
l_rand=(int)(RAND() * LENGTH); /* 交叉位置 */
for (l = 0; l < l_rand; l ++) {
gene[k+SIZE][l] = gene[k][l];
}
for (l = l_rand; l < LENGTH; l ++) {
gene[k+SIZE][l] = gene[k+1][l];
}
for (l = 0; l < l_rand; l ++) {
gene[k+1+SIZE][l] = gene[k+1][l];
}
for (l = l_rand; l < LENGTH; l ++) {
gene[k+1+SIZE][l] = gene[k][l];
}
}
/* 突然変異 */
for (k = SIZE ; k < 2 * SIZE; k ++) {
for (l = 0; l < LENGTH; l ++) {
if (RAND() < PROB) {
gene[k][l] = 1 - gene[k][l];
}
}
}
/* 目的関数の計算 */
for (k = 0; k < 2 * SIZE; k ++) {
total_weight[k] = 0.0;
total_value[k] = 0.0;
for (l = 0; l < LENGTH; l ++) {
total_weight[k] = total_weight[k] + gene[k][l] * weight[l];
total_value[k] = total_value[k] + gene[k][l] * value[l];
}
if (total_weight[k] > LIMIT) {
total_value[k] = 0.0;
}
}
/* ソート */
for (k = 0; k < 2 * SIZE; k ++) {
v_sort[k] = k;
}
for (k = 0; k < 2 * SIZE; k ++) {
for (k_sort = k + 1; k_sort < 2 * SIZE; k_sort ++) {
if (total_value[k] < total_value[k_sort]) {
swap = total_value[k_sort];
total_value[k_sort] = total_value[k];
total_value[k] = swap;
i_swap = v_sort[k_sort];
v_sort[k_sort] = v_sort[k];
v_sort[k] = i_swap;
}
}
}
/* 新世代におきかえ */
for (k = 0; k < 2 * SIZE; k ++) {
for (l = 0 ;l < LENGTH; l ++) {
next_gene[k][l] = gene[v_sort[k]][l];
}
}
for (k = 0; k < 2 * SIZE; k ++) {
for (l = 0; l < LENGTH; l ++) {
gene[k][l] = next_gene[k][l];
}
}
/* 表示 */
gotoxy(1,1); printf(" %5d ",g);
for (l = 0; l < LENGTH; l ++) {
printf("%1d", gene[k][l]);
}
printf(" %8.3f %8.3f\n", total_weight[v_sort[0]], total_value[0]);
LineTo(0, g, total_weight[v_sort[0]]);
LineTo(1, g, total_value[0]);
}
getch();
closegraph();
return 0;
}
/* s1050.c
* ナップザック問題
* (C) H.Ishikawa 1994 2018
*/
#include <stdio.h>
#include "compati.h"
#include "window.h"
#define LENGTH 40 /* 染色体の長さ */
#define SIZE 10 /* 集団の大きさ */
#define MAX_WEIGHT 10.0 /* 品物の重量の最大値 */
#define MAX_VALUE 10.0 /* 品物の価値の最大値 */
#define LIMIT 100.0 /* ナップザックに入る重量の限界値 */
#define G_END 2000 /* 繰り返し数 */
#define PROB 0.01 /* 突然変異の確率 */
#define PROB2 0.5 /* 初期化の確率 */
int gene[2*SIZE][LENGTH] = {0}; /* 染色体 */
int next_gene[2*SIZE][LENGTH] = {0}; /* 新世代の染色体 */
double weight[LENGTH]; /* 品物の重量 */
double value[LENGTH]; /* 品物の価値 */
double total_weight[2*SIZE]; /* 重量の合計 */
double total_value[2*SIZE]; /* 価値の合計 */
int v_sort[2*SIZE]; /* 順番 */
long g; /* 世代を表すforのカウンタ */
int k,l; /* gene[k][l]のためのforのカウンタ */
int l_rand; /* 交叉位置 */
int k_sort; /* ソート用のforのカウンタ */
int i_swap; /* スワップ用変数 */
double swap; /* スワップ用変数 */
int main(void)
{
/* グラフィックの準備 */
opengraph();
SetWindow(0, 0.0,0.0,G_END,120,
SPACE,HIGHT/2-SPACE/2,WIDTH-SPACE,SPACE);
Axis(0, "generation", G_END / 10, "weight", 20);
MoveTo(0, 0.0, 0.0);
SetWindow(1, 0.0,0.0,G_END,250,
SPACE,HIGHT-SPACE,WIDTH-SPACE,HIGHT/2+SPACE/2);
Axis(1, "generation", G_END / 10, "value", 50);
MoveTo(1, 0.0, 0.0);
setlinestyle(0, 0, 2);
srand(3);
/* 染色体のイニシャライズ */
for (k = 0 ; k < SIZE; k ++) {
for (l = 0; l < LENGTH; l ++) {
if (RAND() < PROB2){
gene[k][l] = 1;
}
}
}
/* 品物の重さと価値の設定 */
for (l = 0; l < LENGTH; l ++) {
weight[l] = MAX_WEIGHT * RAND();
value[l] = MAX_VALUE * RAND();
}
/* メイン始まり */
for (g = 1; g <= G_END; g ++) {
/* 交叉 */
for (k = 0; k < SIZE; k = k + 2) {
l_rand=(int)(RAND() * LENGTH); /* 交叉位置 */
for (l = 0; l < l_rand; l ++) {
gene[k+SIZE][l] = gene[k][l];
}
for (l = l_rand; l < LENGTH; l ++) {
gene[k+SIZE][l] = gene[k+1][l];
}
for (l = 0; l < l_rand; l ++) {
gene[k+1+SIZE][l] = gene[k+1][l];
}
for (l = l_rand; l < LENGTH; l ++) {
gene[k+1+SIZE][l] = gene[k][l];
}
}
/* 突然変異 */
for (k = SIZE ; k < 2 * SIZE; k ++) {
for (l = 0; l < LENGTH; l ++) {
if (RAND() < PROB) {
gene[k][l] = 1 - gene[k][l];
}
}
}
/* 目的関数の計算 */
for (k = 0; k < 2 * SIZE; k ++) {
total_weight[k] = 0.0;
total_value[k] = 0.0;
for (l = 0; l < LENGTH; l ++) {
total_weight[k] = total_weight[k] + gene[k][l] * weight[l];
total_value[k] = total_value[k] + gene[k][l] * value[l];
}
if (total_weight[k] > LIMIT) {
total_value[k] = 0.0;
}
}
/* ソート */
for (k = 0; k < 2 * SIZE; k ++) {
v_sort[k] = k;
}
for (k = 0; k < 2 * SIZE; k ++) {
for (k_sort = k + 1; k_sort < 2 * SIZE; k_sort ++) {
if (total_value[k] < total_value[k_sort]) {
swap = total_value[k_sort];
total_value[k_sort] = total_value[k];
total_value[k] = swap;
i_swap = v_sort[k_sort];
v_sort[k_sort] = v_sort[k];
v_sort[k] = i_swap;
}
}
}
/* 新世代におきかえ */
for (k = 0; k < 2 * SIZE; k ++) {
for (l = 0 ;l < LENGTH; l ++) {
next_gene[k][l] = gene[v_sort[k]][l];
}
}
for (k = 0; k < 2 * SIZE; k ++) {
for (l = 0; l < LENGTH; l ++) {
gene[k][l] = next_gene[k][l];
}
}
/* 表示 */
gotoxy(1,1); printf(" %5d ",g);
for (l = 0; l < LENGTH; l ++) {
printf("%1d", gene[k][l]);
}
printf(" %8.3f %8.3f\n", total_weight[v_sort[0]], total_value[0]);
LineTo(0, g, total_weight[v_sort[0]]);
LineTo(1, g, total_value[0]);
}
getch();
closegraph();
return 0;
}
|
１．１０．１節のクラップスは胴元が有利なゲームです。勝っても負けても、毎回１ドルずつ賭け続けた場合、どのくらいの回数でプレーヤが破産するか確かめなさい。
（解答）
| a1010.c １０章１．解答 |
/* a1010.c １０章１．解答
* クラップス
* (C) H.Ishikawa 1994 2018
*/
#include <time.h> /* time() */
#include "compati.h"
#include "window.h"
#define T_END 10000 /* 長期間やってみる */
long t; /* クロック */
long q = 100; /* 持ち金 */
long k = 1; /* １回の賭金 */
int dice(void);
int craps(void);
int main(void)
{
opengraph();
SetWindow(0, 0.0,0.0,(double)T_END,200.0,
SPACE,HIGHT-SPACE,WIDTH-SPACE,SPACE);
Axis(0, "t", T_END/10, "q", 100);
MoveTo(0, 0.0, q);
srand((unsigned)time(NULL)); /* シミュレーションの都度乱数系列を変える */
for (t = 1; t <= T_END; t ++) {
q = q - k;
if (craps() == 1) { /* 勝った場合 */
q = q + 2 * k;
k = 1; /* １ドル賭ける */
}
LineTo(0, (double)(t), (double)(q));
if (q == 0) { break; } /* 破産したら止まる */
}
getch();
closegraph();
return 0;
}
int dice(void)
{
return((int)(RAND() * 6.0 + 1));
}
int craps(void)
{
int r; /* 1:勝ち 0:負け */
int sum1; /* 和１ */
int sum2; /* 和２ */
sum1 = dice() + dice();
if (sum1 == 7 || sum1 == 11) { /* ナチュラル */
r = 1;
} else if(sum1 == 2 || sum1 == 3 || sum1 == 12) { /* クラップス */
r = 0;
} else {
sum2 = 0;
r = 0;
while (sum2 != 7) { /* ７が出るまで */
sum2 = dice() + dice();
if (sum2 == sum1) { /* ポイント */
r = 1;
break;
}
}
}
return(r);
}
|
２． １０．２節の在庫管理問題では、発注点をわった時に発注していたが、毎週１回定時に入荷があるとして、その入荷量を前々日の在庫量により決定し、それにより発注する方法について、シミュレーションし、最適入荷量 を求めなさい。
（解答）
| a1020.c １０章２．解答 在庫管理の問題 （定期入荷） |
/* a1020.c １０章２．解答
* 在庫管理の問題 （定期入荷）
* (C) H.Ishikawa 1994 2018
*/
#include <stdio.h>
#include <math.h>
#include "compati.h"
#include "window.h"
#define FRIDAY 4 /* 入荷日は金曜日 */
#define AL 5 /* １日の平均売上個数 */
#define SP 7000.0 /* 売値 */
#define PC 5000.0 /* 仕入れ単価 */
#define C1 20.0 /* １日１個当たりの在庫費用 */
#define C2 10000.0 /* １回の発注にともなう費用 */
#define C3 500.0 /* 品切れ１個あたりの損失 */
#define T_END 10000
#define A1_MIN 20
#define A1_MAX 200
long poisson(double lambda);
int main(void)
{
long tt; /* シミュレーション日数 */
long a; /* １回の入荷量 */
long q; /* 現在在庫量 */
long a1; /* 発注判断量 */
long k; /* １日の売上個数 */
long zk; /* kの累計 */
long yk; /* 売上個数の累計 */
long ls; /* 在庫割れのため売り損なった個数 */
long nn; /* 発注回数 */
double sc1; /* 在庫費用の累計 */
double sc2; /* 発注費用の累計 */
double sc3; /* 品切れ損失の累計 */
double rr; /* 総売上 */
/* グラフィック準備 */
opengraph();
SetWindow(0, 0.0,0.0,A1_MAX,1.0e08,
SPACE,HIGHT-SPACE,WIDTH-SPACE,SPACE);
Axis(0, "q1", 10, "rr", 1e07);
setlinestyle(0, 0, 2);
printf(" a1 yk ls nn sc1 sc2 sc3 rr\n");
/* 開始 */
for (a1 = A1_MIN; a1 <= A1_MAX; a1 = a1 + 10) {
/* 初期設定 */
srand(113);
tt = 0; q = a1; zk = 0; ls = 0; nn = 0;
sc1 = 0.0; sc2 = 0.0; sc3 = 0.0;
/* １日分 */
while (tt <= T_END) {
k = poisson(AL); /* k個売れた */
q = q - k;
zk = zk + k;
if (q < 0) { /* 在庫割れのため売り損ない */
ls = ls - q;
q = 0;
}
if ((tt % 7) == (FRIDAY - 2)) { /* 入荷量決める日は金曜の２日前 */
a = a1 - q; /* a1と今日の在庫の差が入荷 */
if (a <= 0) {
a = 0;
} else {
nn = nn + 1;
}
}
if ((tt % 7) == FRIDAY) { /* 入荷日 */
q = q + a;
}
sc1 = q * C1 + sc1;
tt = tt + 1;
}
sc2 = nn * C2;
sc3 = ls * C3;
yk = zk - ls;
rr = (SP - PC) * yk - sc1 - sc2 - sc3;
printf("%4ld %6ld %6ld %6ld %8.0f %8.0f %8.0f %8.0f\n",
a1, yk, ls, nn, sc1, sc2, sc3, rr);
if (a1 == A1_MIN) {
MoveTo(0, a1, rr);
} else {
LineTo(0, a1, rr);
}
}
getch();
closegraph();
return 0;
}
long poisson(double lambda)
{
double xp = RAND();
long k = 0;
while (xp >= exp(-lambda)) {
xp = xp * RAND();
k = k + 1;
}
return (k);
}
|
３．プログラムs1040.cで w を wc を中心に変えて実行してみなさい。 w が wc より小さい場合、磁化率が０とならず磁石のままであることをたしかめなさい。また
w が wc より大きい場合、磁化率が０となる速度が早いことをたしかめなさい。
（解答省略）
|
 |
