本文へスキップ
モンテカルロシミュレーション
10. 何でもシミュレーションしてみよう
 本書の最後は、「何でもシミュレーションしてみよう」と題し、ギャンブルから、最近注目されている、ジェネティックアルゴリズムまで、乱数を用いた、さまざまな、確率実験を、試みます。

クラップスという名のゲーム文献(8)(13)
 カジノで行われているダイス(サイコロ)を使った、クラップスというゲームをとりあげます。このクラップスに限らず、カジノで行われるゲームは、すべて 胴元が有利となるようなルールが定められています。したがって長時間平均すれば、必ず胴元が勝つ仕組みとなっているわけですが、必勝法はないものでしょう か。
 まず、クラップスのルールを説明します。クラップスは2個のサイコロを用いておこないます。2個のサイコロの目の合計は2から12までの値をとりま す。(3.5節で示したとおり、7の出る確率が最も大きい。)プレーヤは2個のサイコロを振って、その和が7または、11の時(これをナチュラルといいま す)は勝ちで,2、3あるいは12の時(これをクラップスといいます)は負けです。もしそれ以外のとき、すなわち、4、5、6、8、9または、10がでた ときは、その値をポイントといい、ひきつづき同じポイントまたは7が出るまで、サイコロを振ります。そして同じポイントがでれば、プレーヤの勝ち、それ以 前に7が出れば、負けとなります。もしプレーヤが勝てば、自分の賭け金の2倍の払い戻しを受けます。
 このゲームは、たいへんよくできて、ごく僅かな確率の差で、胴元が勝つようにできています。


必勝法
 このゲームで必ずプレーヤが勝つ方法はあるでしょうか。もしそれがあれば、ぜひカジノに行って見たいものです。倍増法という必勝法が知られています。倍 増法というのは、まず1ドルの賭け金からスタートし、もし勝てば、次も1ドルを賭けますが、負けた場合は2ドルを、さらに負ければ前回の賭け金の2倍の額 を次々に賭けていく方法です。この方法では、負けた金額をとりもどせるので、必勝法となりそうです。はたして、これがほんとうに必勝法となるのでしょう か。クラップスにおいて、倍増法をシミュレーションしてみましょう。

シミュレーション・プログラム
 ではクラップスのシミュレーションプログラムs1010.cの説明を行ないます。まず、サイコロの関数を作ります。
   int dice(void);

により、整数1〜6の一様分布の乱数を発生させます。次に、(2)のとおりクラップスのルールどおりに記述して、関数
   int craps(void);

を作ります。1回目のサイコロの目の和を sum1 , 2回目以後のサイコロの目の和を sum2 としています。この関数は、勝ちの場合1、負けの場合0の値を返します。
 次にメインプログラムですが、1回の賭け金を k ドル、現在の持ち金を q ドルとします。行30から行35が倍増法による賭けで、もし勝った場合、次の賭け金を1とし、負けた場合前回の2倍を賭けることとします。プログラムでは、賭けを何回も繰り返し、横軸賭けの回数 t 、縦軸持ち金 q をグラフに表示するようにしています。

シミュレーションの結果
 プログラムを実行すると、図10.1.1のように時間とともに、着実に持ち金は増加していきますが、時々負けが続くと賭け金が急増し、赤字に転落すると きがあります。もしその時、もち金がなくなれば破産ですが、なお余裕があって倍々と賭けることができれば、必ずプラスになります。しかしながら、通常は胴 元のほうが資金は豊富で、負けが続いても持ちこたえられるのに対し、客のほうが先に底をついてしまいます。もし十分な資金を用意することができれば、ぜひ カジノへ行って、ひと儲けして下さい。

図10.1.1 倍増法による必勝法

  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;                  /* 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;            /* 1ドル賭ける */
    } else {            /* 負けた場合 */
      k = k * 2;        /* 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;         /* 和1 */
  int     sum2;         /* 和2 */
  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) {             /* 7が出るまで */
      sum2 = dice() + dice();
      if (sum2 == sum1) {           /* ポイント */
        r = 1;
        break;
      }
    }  
  }
  return(r);    
}


在庫管理問題とは
 オペレーションズリサーチの教科書に関らず出てくる問題です。企業や産業システムでは、原料や製品が供給から需要ひ向かって流れていく過程において、それらを貯蔵する機能が必要です。この貯蔵容量と、それを運用管理する最適の方法を決定するのが、在庫管理問題です。
 例えば、ある商店が商品を問屋から仕入れて、客に販売する場合を考えてみましょう。商品を一度に多量仕入れると、在庫をたくさんかかえ、そのための費用 がかさみます。逆に少量の在庫では、大きな需要があった場合に品切れの危険が生じ、得べかりし利益をそこない、ひいてはその店の信用を失います。しから ば、@ 在庫に要する費用、A 問屋への発注のための費用、B 品切れが発生した場合にこうむる損失の合計を最小にして、利益を最大にするには、いつどれ だけ問屋に発注すべきでしょうか。
 販売店の在庫 q は、図10.2.1のように毎日の売上げ k 個ずつ減少していきます。基準発注点 q1 をわったら、 a 個の発注をしますが、すぐには入荷せず、 lt 日の遅れがあります。図10.2.2のように、在庫量は変化するでしょう。売上げ個数と入荷おくれが確率的である場合を考え、 基準発注点q1 の最適値を求めるためのシミュレーション・プログラムを作ってみましょう。

図10.2.1 在庫モデル

図10.2.2 在庫 q の変化と lt



プログラムのポイント
 売上げ個数は、毎日平均 AL のポアソン分布に従うとします。また入荷遅れ日数 lt は平均 ET 、標準偏差 ST の正規分布にしたがうものとします。これらの乱数の発生には4章で作った関数を用います。シミュレーション・プログラムの方式は1日をクロックとして進める固定時間方式を用います。プログラムs1020.cの行60から行78が,1日分を疑似します。ここで重要な変数は、入荷遅れ日数 lt で次の意味を持ちます。


   lt > 0 :すでに発注ずみであり、あと lt 日で納品される。
   lt = 0 :本日納品される。
   lt < 0 :発注していない。
 行61で1日平均 k 個ポアソン分布の売上げがあったとします。行64で発注点 q1 を割り、まだ発注していなければ発注します。行68で q が負になれば、在庫なしで販売しそこなった数 ls を数えます。行73では、本日納品されたので、在庫に A を加えます。行75で在庫費用の累計を計算します。 T_END 日分のシミュレーションを終わったあと、発注費用の累計 sc2 、品切れ損失の累計 sc3 、売上個数の累計 yk を計算し、利益 rr を求めます。そしてグラフをプロットし、発注点 q1 をかえ、シミュレーションを続けます。

アウトプット
 10000日分のシミュレーションを q1 を0から100まで変化させ、繰り返し行います。 発注点 q1に対する総売上げ rrのグラフ図10.2.3を見ると、q1=45 の時に、利益が最大となることがわかります。

図10.2.3 在庫管理の問題 シミュレーション結果

  s1020.c   在庫管理の問題
/* s1020.c 
 *        在庫管理の問題 
 *         (C) H.Ishikawa  1994  2018
 */

#include <stdio.h>
#include <math.h>
#include "compati.h"
#include "window.h"

#define A        100                /* 1回の入荷量 */
#define AL       5                  /* 1日の平均売上個数 */
#define ET       6                  /* 納品するまでの平均日数 */
#define ST       2                  /* 納品するまでの日数の標準偏差 */
#define SP       7000.0             /* 売値 */
#define PC       5000.0             /* 仕入れ単価 */
#define C1       20.0               /* 1日1個当たりの在庫費用 */
#define C2       10000.0            /* 1回の発注にともなう費用 */
#define C3       500.0              /* 品切れ1個あたりの損失 */
#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;                   /* 1日の売上個数 */
  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;
    
    /* 1日分のシミュレーション */
    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);
}


鉄は強磁性体
 本節では、シミュレーション技術を用いて、1つ1つの原子の動きを計算し、物質全体としてどういう性質をもっているかを、鉄のような強磁性体を例に解析してみましょう。
 原子の模型は、図10.4.1のように、原子核のまわりを電子がまわっているように書けます。電子はマイナスの電気をもっており、電子が回転運動すると 磁場ができます。その軸は原子の中心を通り、電子の軌道面に対して垂直です。この磁場のことを電子の軌道磁気モーメントと言います。また電子は自転しなが ら軸のまわりを回っており、これによっても、磁場が生じています。これをスピン磁気モーメントと言います。これら二種類の磁気モーメントが加え合わされ、 原子1つの磁気モーメントとなります。
図10.4.1 原子の模型


 鉄の原子の場合は、スピン磁気モーメントの釣合がひどくずれており、原子1つ1つが小さな磁石となっており、その強さはすべての元素の原子なかで一番強 いのです。鉄の原子は、結晶の格子点に固定されて動くことができませんが、その向きは自由に変り、いろいろの方向を向きます。鉄の棒を磁化して磁石にする というのは、できるだけたくさんの原子の磁場の方向を並行にそろえることにほかなりません。
 これはごく大ざっぱにいうと、鉄原子の電子スピンによる磁気モーメントは、近接する原子間では、互いに逆向きよりも平行の方がエネルギーが低くなるような相互作用が動いているためです。(図10.4.2a)


磁石が磁石でなくなる
 ところが温度を高くすると、原子はしだいに無秩序になり、上向き下向きのスピンの原子があらわれます。そしてある温度(ピエル・キュリーが発見した。 キュリー温度という)を越すと、向きがばらばらになり、鉄全体の磁気モーメントはほぼ零になり、鉄は磁石でなくなってしまうのです(図10.4.2b)。 このことを、強磁性体からの相転移といいます。シミュレーションにより、この相転移の実験をしてみましょう。
図10.4.2 絶対零度の場合とある温度の場合のスピン磁気モーメント


イジング・モデル
 今原子が図10.4.3のように2次元格子状に並んでいるとします。原子 (i,j) の状態は、スピンが上向きか、下向きかという2つの状態だけをとるものとし、それを s(i, j)=1 , s(i, j)= -1 とします。この値をスピン変数といいます。
 ある原子 (i0, j0) が、まわりの原子から受けるエネルギーは、まわりの原子のスピンの方向と、自分のスピンの方向によって決まり、その値は、
・・・・・ 式10.4.1
となります。 J は物質によって異なり、 (i0, j0) から離れるにしたがって小さな値となります。特に原子 (i0, j0) から見てすぐ隣りの原子からのみ相互作用を受けるとすれば、式は、
 ・・・式10.4.2    
となります。 e の値は、原子 (i0, j0) から見て、まわりの原子のスピンの方向が全部逆であれば、一番大きく、また原子 (i0, j0) とまわりの原子のスピンが全部同じであれば、一番小さくなります。
 ある温度で、原子 (i0, j0) が向きを反転する確率は、統計力学の教えるところにより、単位時間あたり、
・・・・・ 式10.4.3
となります。ここで、 k はボルツマン定数、 T は絶対温度、 e は式の値、 A は定数です。
 式10.4.2と、式10.4.3から、原子 (i0, j0) が反転する確率は、 T が高くなると高くなり、逆向きの原子にとりかこまれているとき最も高く、まわりの原子と同じ向きのとき最も低くなります。
 このようなモデルを提案者の名前をとって、イジング・モデルといいます。

図10.4.3


プログラムのポイント
 以上のモデルを、シミュレーションにより実現します(プログラムs1040.c)。温度を高温にすると、あちらこちらで、同時に反転が起こりますが、シミュレーションでは、時計を停止させ、ランダムに原子を選び出し、処理します。
 それにはまず、モデルの磁石の原子に相当する配列を用意します。行14で1辺の大きさが、 MAX+2 (値は 0 から MAX+1 )の配列を作り、その中の 1 から MAX までを使います。 0MAX+1 は式10.4.2の計算をするとき ij1MAX の値をとったとき、その隣りの値が配列をはみださないようにするためです。そして初期状態は、絶対零度としすべて磁気モーメントは同方向、すなわち s(i,j)=1 とします。
 次に行47,行48で2つ1組のから MAX までの整数値をとる一様分布の乱数を発生させます。この1組の乱数により、原子を1つ選び、式10.4.2の相互作用のエネルギーを計算し(行49)、式10.4.3の反転の確率を計算します(行50)。ここで w0 = kT/J とおいて正規化温度としています。また A = exp(-4/w0) とします。別に発生した乱数とこの確率を比べ(行51)、乱数の方が小さければ、原子を反転させます(行52)。この処理を MAX×MAX 回繰り返すと、各原子について、平均1回処理を行ったことになるので、クロックを1進めます。

アウトプット
  s(i, j)=1 である原子の数を数えて、磁化率を計算します(行62)。磁化率は全部の原子が s(i, j)=1 であれば1, 1と−1が半々であれば0の値を取ります。また、ディスプレイに、刻々と原子の反転の様子をプロットさせます。それには反転のたびに s(i, j) が1の時は青、−1の時は赤で表示します。これによりミクロな原子の動きをパターンとして視覚に訴えることができます。
 プログラムを実行してみましょう。図10.4.4のように絶対零度では、一面に正のモーメント(図では青)であったものが、しだいに負のモーメントの部分(図では赤)ができ、時間の経過とともにそれが合体し(これを磁気ドメインといいます)、形も刻々変化していきます。
 ウィンドウ1のグラフは、磁化率をあらわします。
 ある温度(キュリー点)以下では、磁化率は0にならず磁石のままで、ある温度 Tc をこえると、相転移が起こり磁化率が0となり、磁石が磁石でなくなっていきます。この時の正規化温度の、理論値は、
・・・・・ 式10.4.4
で与えられます。 プログラムでは正規化温度w0 を変えることができるようになっていますから、 wc を中心に色々変えて相転移が起こる様子を実験してみて下さい。


図10.4.4 イジングモデル シミュレーション結果
しだいに負のモーメントの部分が増加(図では赤)
時間の経過とともに、うごめうごめきながらそれが合体していく


  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;
}




ジェネティック・アルゴリズムとは
 ジェネティック・アルゴリズム(GA)は遺伝的アルゴリズムともよばれ、生物の進化過程をシミュレーションすることによって、たくさんの組み合せのなか から最適な解を見つけ出す手法のことをいいます。ダーウィンの自然淘汰理論というのがありますが、これを基礎に、生物の進化過程を工学的モデルにおきかえ たものです。
 進化過程の概要は、おおよそ次のようにあらわされるでしょう。まず、0と1の値をもつ遺伝子(Gene)からなる染色体(Chromosome)を考え ます。この染色体のあつまりを集団(Population)といいます。それに対し図10.5.1ように(1)、(2)、(3)の操作をほどこします。
(1) 交叉
 集団内の2つの染色体の組(ペア)において、交叉(crossover)という遺伝子の一部の入れ替えがおこなわれること。
(2) 突然変異
 各染色体において、ある確率で、遺伝子の一部が0から1あるいは、1から0に返転すること。
(3) 選択
 (1)、(2)をへた遺伝子をもった染色体のうち環境へのより高い適応性をもつものが、選択されること。
(1)、(2)、(3)の遺伝子的操作を1世代とよび、この操作をくりかえしていくと、その集団はより環境に順応する方向に進化します。
 ジェネティック・アルゴリズムでは、さまざまな組み合わせを遺伝子からなる染色体で表現し、その組み合わせごとの目的関数を計算します。そして、その値 の高いものを選択し、それを残すよう生存競争を何世代にもわたり行わせれば、最終的には、目的関数の値が最大となるようにすることができます。
 ジェネティック・アルゴリズムの原理は、驚くほど単純ですが、その探索能力は非常に高く、通常の方法では解を見つけにくい場合でも、比較的容易に正確に近づくことができるといわれています。

図10.5.1 組み替え(交叉と突然変異)


ナップザック問題
 ジェネティック・アルゴリズムのメカニズムを「ナップザック問題」という典型的な組み合せ問題を例に、プログラムを作ってみましょう。
 図10.5.2のように、ナップザックがあります。ナップザックにはこべる重さには限界があります。たくさんの異なる重さの、異なる価値をもつ品物が与 えられたとき、重さがナップザックの限界値内で品物を幾つか選択し、その時の価値の合計が最大となるような組み合せを見つけるという問題です。すなわち、 この場合の目的関数は価値の合計で、ナップザックがこわれない範囲で、価値が最も高くなるよう、品物を選んで運ぼうとするものです。通常の方法では解を見つけにくい例です。
図10.5.2 ナップザックの問題


ナップザック問題のプログラム
 ナップザック問題を、ジェネティック・アルゴリズムの手法で解く場合、まず、ナップザックに品物が入っている状態を染色体に対応させて、表現する必要があります。図10.2.5の例では、品物1、3、4をナップザックに入れている状態ですから、これを染色体で
    01011
と表現します。
 プログラムs1050.cについて説明していきます。LENGTHという長さをもつ染色体をSIZE個用意します。実際には組み替え操作のため,この2倍の配列を行19により用意します。また,行20により新世代の染色体を表わす配列も準備します。このほか1つ1つの品物の重量,価値をあらわす配列,染色体に対応した重量の合計,価値の合計をあらわす配列を用います。
 行50から行56までで染色体を初期状態にします。PROB2の確率で1にセットし,また行59から行62において,各品物ごとの重量と価値も乱数でセットしておきます。
 行65からがメインの始まりで,遺伝子の組み替え(交叉と突然変異),目的関数の計算,ソート,新世代へのおきかえ,表示という手続きを,G_ENDで示される世代にわたって繰り返し計算します。
 遺伝子の組み替えは,図10.5.3のように行います。k番目の染色体と,k+SIZE番目の染色体がペアであるとして,交叉の操作を行います。行69のように交叉位置を乱数で発生させ,k番目の後半とk+SIZE番目の後半を入れ換えます。次に行85以下のようにkの値がSIZE以上の染色体についてPROBという確率で,1と0を反転させ,突然変異を起こさせます。
 次に,目的関数を計算します。目的関数は,染色体ごとの,合計の重量と合計の価値です。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]を入れ換えます。
 
図10.5.3 遺伝子の組み換え



結果の出力
 2つのウインドウに横軸世代、縦軸それぞれ重量の合計、価値の合計をグラフ表示します。また、もっとも目的関数の大きい k = 0 の染色体の構造、またそのときの重量の合計、価値の合計を表示します。
 このプログラムを走らせると、図10.5.4のように、100世代ほどで、最適に近い値が得られていますが、1000世代くりかえし実行する必要があります。

図10.5.4 ナップザック問題 シミュレーション結果

  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;
}


.10.1節のクラップスは胴元が有利なゲームです。勝っても負けても、毎回1ドルずつ賭け続けた場合、どのくらいの回数でプレーヤが破産するか確かめなさい。
  (解答) 
  a1010.c 10章1.解答
/* a1010.c    10章1.解答
 *        クラップス
 *         (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;                  /* 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;                    /* 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;         /* 和1 */
  int     sum2;         /* 和2 */
  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) {             /* 7が出るまで */
      sum2 = dice() + dice();
      if (sum2 == sum1) {           /* ポイント */
        r = 1;
        break;
      }
    }  
  }
  return(r);    
}

. 10.2節の在庫管理問題では、発注点をわった時に発注していたが、毎週1回定時に入荷があるとして、その入荷量を前々日の在庫量により決定し、それにより発注する方法について、シミュレーションし、最適入荷量 を求めなさい。
  (解答) 
  a1020.c 10章2.解答   在庫管理の問題  (定期入荷)
/* a1020.c   10章2.解答
 *        在庫管理の問題  (定期入荷)
 *         (C) H.Ishikawa  1994 2018
 */

#include <stdio.h>
#include <math.h>
#include "compati.h"
#include "window.h"

#define FRIDAY   4                  /* 入荷日は金曜日 */
#define AL       5                  /* 1日の平均売上個数 */
#define SP       7000.0             /* 売値 */
#define PC       5000.0             /* 仕入れ単価 */
#define C1       20.0               /* 1日1個当たりの在庫費用 */
#define C2       10000.0            /* 1回の発注にともなう費用 */
#define C3       500.0              /* 品切れ1個あたりの損失 */
#define T_END    10000
#define A1_MIN    20
#define A1_MAX    200

long poisson(double lambda);

int main(void)
{
  long    tt;                  /* シミュレーション日数 */
  long    a;                   /* 1回の入荷量 */
  long    q;                   /* 現在在庫量 */
  long    a1;                  /* 発注判断量 */
  long    k;                   /* 1日の売上個数 */
  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;
    
    /* 1日分 */
    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)) {   /* 入荷量決める日は金曜の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 より小さい場合、磁化率が0とならず磁石のままであることをたしかめなさい。また w が wc より大きい場合、磁化率が0となる速度が早いことをたしかめなさい。
  (解答省略)


top