GA により 2 変数関数の最小値を求める

データ例，プログラム（ MT.h は除く）

　　以下，複数のファイル構成になっています．ファイル間の区切りを「---･･･」で示します．このプログラム例においては，ヘッダファイル MT.h に記述されたメルセンヌ・ツイスタを使用していまが，C++11 で記述可能であれば，C++ の標準ライブラリであるメルセンヌ・ツイスター法による擬似乱数生成エンジンを使用できます． 

-------------------------i_data-------------------------
世代交代数 1000 乱数の初期値 12345 関数 3 出力 100
個体数 20 子供の数 20 遺伝子長 20
エリート 2 突然変異率 0.03

-------------------------program-------------------------
/*************************/
/* クラス Species の定義 */
/*************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "MT.h"

class Species {
		double *pi;   // 適応度
		double *ro;   // ルーレット板
		int size;   // 個体総数
		int max_ch;   // 子供の数の最大値
		int len;   // 遺伝子長
		unsigned char **ind;   // 集団（個体の集まり）
		char *pi_w;   // 適応度計算指標
                      //   =0 : 未使用
                      //   =1 : 適応度計算前（突然変異はこの個体だけに適用）
                      //   =2 : 適応度計算済み（交叉時に親とみなす）
		char *s_w;   // 淘汰用指標（選択された回数）
	public:
                    // コンストラクタ
		Species(int, int, int, long);
                    // デストラクタ
		~Species();
                    // 出力
		void O_std(int, int out=20);
                    // 初期設定
		void I_std();
                    // 交叉
		void C_copy();   // 親のコピー
		void C_unifm(int pair=0, double pr=1.0);   // 一様交叉
                    // 突然変異
		void M_alle(double);   // 対立遺伝子への置換
                    // 淘汰
		void S_roul(int elite=2);
                    // 評価（最大適応度を返す）
		double Adap(int);   // 個体の適応度の計算
                    // その他
		int position();   // 場所を探す
		int select();   // 親の選択
};

/****************************/
/* コンストラクタ           */
/*      s : 個体総数        */
/*      mc : 子供の数       */
/*      mx : 遺伝子長       */
/*      seed : 乱数の初期値 */
/****************************/

Species::Species(int s, int mc, int mx, long seed)
{
	int i1, num;
/*
     乱数の初期設定
*/
	init_genrand(seed);
/*
     値の設定
*/
	size   = s;
	max_ch = mc;
	len    = mx;
/*
     領域の確保
*/
	num = size + max_ch;

	ind = new unsigned char * [num];
	for (i1 = 0; i1 < num; i1++)
		ind[i1] = new unsigned char [2*len];

	pi   = new double [num];
	pi_w = new char [num];
	s_w  = new char [num];
	ro   = new double [num];
}

/****************/
/* デストラクタ */
/****************/
Species::~Species()
{
	int i1;

	for (i1 = 0; i1 < size+max_ch; i1++)
		delete [] ind[i1];
	delete [] ind;

	delete [] pi;
	delete [] pi_w;
	delete [] s_w;
	delete [] ro;
}

/****************************/
/* ２変数関数の最小値       */
/*      coded by Y.Suganuma */
/****************************/
int main(void)
{
	double mute;
	long seed;
	int fun, gen, max_gen, o_lvl, size, child, len, elite;
					// データの入力
	scanf("%*s %d %*s %ld %*s %d %*s %d", &max_gen, &seed, &fun, &o_lvl);
	scanf("%*s %d %*s %d %*s %d", &size, &child, &len);
	scanf("%*s %d %*s %lf", &elite, &mute);
					// 種の定義
	Species a(size, child, len, seed);
					// 初期設定
	gen = 0;
	a.I_std();
                    // 個体の評価
	a.Adap(fun);
					// 初期世代の出力
	if (o_lvl > 0)
		a.O_std(gen);
					// 実行
	for (gen = 1; gen <= max_gen; gen++) {
						// 交叉
		a.C_unifm();
						// 突然変異
		a.M_alle(mute);
						// 適応度の計算
		a.Adap(fun);
						// 淘汰
		a.S_roul(elite);
						// 出力
		if (o_lvl > 0 && gen%o_lvl == 0)
			a.O_std(gen);
	}
					// 最終出力
	if (o_lvl == 0)
		a.O_std(gen-1, size+1);

	return 0;
}

/************************/
/* 適応度の計算         */
/*      sw : 関数の種類 */
/************************/
double Species::Adap(int sw)
{
	double cv = 5.0 / (pow(2.0, (double)len) - 1.0), max = 0.0, x, x1, x2, x3, y, w;
	int i1, i2, k = -1;

	for (i1 = 0; i1 < size+max_ch; i1++) {
		if (pi_w[i1] == 1) {

			x = 0.0;
			w = 0.0;
			for (i2 = len-1; i2 >= 0; i2--) {
				if (ind[i1][i2] > 0)
					x += pow(2.0, w);
				w += 1.0;
			}
			x *= cv;

			y = 0.0;
			w = 0.0;
			for (i2 = 2*len-1; i2 >= len; i2--) {
				if (ind[i1][i2] > 0)
					y += pow(2.0, w);
				w += 1.0;
			}
			y *= cv;

			pi_w[i1] = 2;
			switch (sw) {
				case 1:
					x1     = x - 1.0;
					x2     = y - 2.0;
					pi[i1] = -x1 * x1 - x2 * x2;
					break;
				case 2:
					x1     = y - x * x;
					x2     = 1.0 - x;
					pi[i1] = -100.0 * x1 * x1 - x2 * x2;
					break;
				case 3:
					x1     = 1.5 - x * (1.0 - y);
					x2     = 2.25 - x * (1.0 - y * y);
					x3     = 2.625 - x * (1.0 - y * y * y);
					pi[i1] = -x1 * x1 - x2 * x2 - x3 * x3;
					break;
			}
		}

		if (pi_w[i1] > 0 && (k < 0 || pi[i1] > max)) {
			max = pi[i1];
			k   = 1;
		}
	}

	return max;
}

/************/
/* 初期設定 */
/************/
void Species::I_std()
{
	int i1, i2, i3, sw1, sw2;
/*
     適応度
*/
	for (i1 = 0; i1 < size+max_ch; i1++) {
		if (i1 < size)
			pi_w[i1] = 1;   // 適応度の計算前
		else
			pi_w[i1] = 0;   // 未使用
	}
/*
     染色体
*/
	for (i1 = 0; i1 < size; i1++) {
		sw1 = 0;
		while (sw1 == 0) {
                              // 遺伝子の決定
			for (i2 = 0; i2 < 2*len; i2++)
				ind[i1][i2] = (genrand_real3() > 0.5) ? 1 : 0;
                              // 重複個体のチェック
			sw1 = 1;
			for (i2 = 0; i2 < i1 && sw1 > 0; i2++) {
				sw2 = 0;
				for (i3 = 0; i3 < 2*len && sw2 == 0; i3++) {
					if (ind[i1][i3] != ind[i2][i3])
						sw2 = 1;
				}
				if (sw2 == 0)
					sw1 = 0;
			}
		}
	}
}

/*********************************************/
/* 標準出力                                  */
/*      gen : 現在の世代                     */
/*      out : n>0 : ｎ個ずつ出力(default=20) */
/*********************************************/
void Species::O_std(int gen, int out)
{
	double cv = 5.0 / (pow(2.0, (double)len) - 1.0), x, y, w;
	int count = 0, i1, i2, num = 0;

	printf("第 %d 世代\n", gen);

	for (i1 = 0; i1 < size+max_ch; i1++) {
		if (pi_w[i1] > 0) {

			num++;
			printf("%3d ", num);

			x = 0.0;
			w = 0.0;
			for (i2 = len-1; i2 >= 0; i2--) {
				if (ind[i1][i2] > 0)
					x += pow(2.0, w);
				w += 1.0;
			}
			x *= cv;

			y = 0.0;
			w = 0.0;
			for (i2 = 2*len-1; i2 >= len; i2--) {
				if (ind[i1][i2] > 0)
					y += pow(2.0, w);
				w += 1.0;
			}
			y *= cv;

			printf("%8.5f %8.5f", x, y);

			if (pi_w[i1] > 1)
				printf(" 適応度 %8.5f\n", pi[i1]);
			else
				printf("\n");

			count++;
			if (count == out) {
				count = 0;
				getchar();
			}
		}
	}

	if (count > 0)
		getchar();
}

/**************/
/* 親のコピー */
/**************/
void Species::C_copy()
{
	int i1, i2, k, p, p1, p2, sw;
                         // 親の選択
	p1 = select();
	sw = 0;
	while (sw == 0) {
		p2 = select();
		if (p1 != p2)
			sw = 1;
	}
                         // コピー
	for (i1 = 0; i1 < 2; i1++) {
		p       = (i1 == 0) ? p1 : p2;
		k       = position();
		pi_w[k] = 1;
		for (i2 = 0; i2 < 2*len; i2++)
			ind[k][i2] = ind[p][i2];
	}
}

/*******************************************************/
/* 交叉（一様交叉．[0,1]を等確率で発生させ，1であれば，*/
/*       親１，0であれば親２の遺伝子を子１が受け継ぐ） */
/*      pair : 親のペア数(default=最大子供数/2)        */
/*      pr : 交叉確率(default=1.0)                     */
/*******************************************************/
void Species::C_unifm(int pair, double pr)
{
	int i1, i2, k1, k2, p1, p2, sw;
/*
     初期設定
*/
	if (pair == 0)
		pair = max_ch / 2;
/*
     交叉
*/
	for (i1 = 0; i1 < pair; i1++) {
                         // 交叉しない場合
		if (genrand_real3() > pr)
			C_copy();
                         // 交叉する場合
		else {
                              // 親の選択と子供を入れる場所
			p1 = select();
			sw = 0;
			while (sw == 0) {
				p2 = select();
				if (p1 != p2)
					sw = 1;
			}
			k1       = position();
			pi_w[k1] = 1;
			k2       = position();
			pi_w[k2] = 1;
                              // 交叉
			for (i2 = 0; i2 < 2*len; i2++) {
				if (genrand_real3() > 0.5) {
					ind[k1][i2] = ind[p1][i2];
					ind[k2][i2] = ind[p2][i2];
				}
				else {
					ind[k1][i2] = ind[p2][i2];
					ind[k2][i2] = ind[p1][i2];
				}
			}
		}
	}
}

/**************************************/
/* 突然変異（対立遺伝子との置き換え） */
/*      pr : 突然変異率               */
/**************************************/
void Species::M_alle(double pr)
{
	int i1, i2;

	for (i1 = 0; i1 < size+max_ch; i1++) {
		if (pi_w[i1] == 1) {
			for (i2 = 0; i2 < 2*len; i2++) {
				if (genrand_real3() <= pr)
					ind[i1][i2] = (genrand_real3() > 0.5) ? 1 : 0;
			}
		}
	}
}

/************************************************/
/* 淘汰（エリート・ルーレット選択）             */
/*      elite : エリートで残す個体数(default=2) */
/************************************************/
void Species::S_roul(int elite)
{
	int i1, i2, i3, k = 0, max, n = 0, p, sw;
/*
     初期設定
*/
	for (i1 = 0; i1 < size+max_ch; i1++)
		s_w[i1] = 0;
/*
     重複個体を削除
*/
	for (i1 = 0; i1 < size+max_ch; i1++) {
		if (pi_w[i1] > 0) {
			for (i2 = i1+1; i2 < size+max_ch; i2++) {
				if (pi_w[i2] > 0) {
					sw = 0;
					for (i3 = 0; i3 < 2*len && sw == 0; i3++) {
						if (ind[i1][i3] != ind[i2][i3])
							sw = 1;
					}
					if (sw == 0)
						pi_w[i2] = 0;
				}
			}
		}
	}

	for (i1 = 0; i1 < size+max_ch; i1++) {
		if (pi_w[i1] > 1)
			n++;
	}

	if (n == 0) {
		printf("***error  残す個体がない (S_roul)\n");
		exit(1);
	}
/*
     淘汰して残す個体を選ぶ
*/
                         // エリートの選択
	sw = 0;

	while (k < elite && k < n && sw == 0) {
		max = -1;
		for (i1 = 0; i1 < size+max_ch; i1++) {
			if (pi_w[i1] > 1 && s_w[i1] == 0) {
				if (max < 0 || pi[i1] > pi[max])
					max = i1;
			}
		}
		if (max < 0)
			sw = 1;
		else {
			s_w[max] = 1;
			k++;
		}
	}
                         // ルーレット選択
	int count = 0;

	while (count < size+max_ch && k < size && k < n) {
		p = select();
		if (s_w[p] > 0) {
			sw = 1;
			count++;
		}
		else {
			count = 0;
			s_w[p]++;
			k++;
		}
	}
                                   // 選択に失敗した場合の処理
	if (k < size && k < n) {
		for (i1 = 0; i1 < size+max_ch && k < size && k < n; i1++) {
			if (pi_w[i1] > 1 && s_w[i1] == 0) {
				s_w[i1] = 1;
				k++;
			}
		}
	}
                              // 複数回選択されたものの処理
	for (i1 = 0; i1 < size+max_ch; i1++) {
		if (s_w[i1] == 0)
			pi_w[i1] = 0;
	}

	for (i1 = 0; i1 < size+max_ch; i1++) {
		if (s_w[i1] > 0) {
			if (s_w[i1] > 1) {
				for (i2 = 2; i2 <= s_w[i1]; i2++) {
					k       = position();
					pi_w[k] = 2;
					pi[k]   = pi[i1];
					for (i3 = 0; i3 < 2*len; i3++)
						ind[k][i3] = ind[i1][i3];
				}
			}
		}
	}
}

/************************/
/* 空いている場所を探す */
/************************/
int Species::position()
{
	int i1, k = -1;

	for (i1 = 0; i1 < size+max_ch && k < 0; i1++) {
		if (pi_w[i1] == 0)
			k = i1;
	}

	if (k < 0) {
		printf("***error  空いている場所がない --position--\n");
		exit(1);
	}

	return k;
}

/**************/
/* 個体の選択 */
/**************/
int Species::select()
{
	double sum, x;
	int i1, k, n = 0, sw;
/*
     ルーレット版の用意
*/
	for (i1 = 0; i1 < size+max_ch; i1++) {
		if (pi_w[i1] > 1)
			n++;
	}
	sum = 1.0 / n;
	for (i1 = 0; i1 < size+max_ch; i1++) {
		if (pi_w[i1] > 1)
			ro[i1] = sum;
	}

	sum = 0.0;
	for (i1 = 0; i1 < size+max_ch; i1++) {
		if (pi_w[i1] > 1) {
			sum    += ro[i1];
			ro[i1]  = sum;
		}
	}
/*
     選択
*/
	x  = genrand_real3();
	sw = 0;
	k  = 0;
	for (i1 = 0; i1 < size+max_ch && sw == 0; i1++) {
		if (pi_w[i1] > 1) {
			if (x <= ro[i1]) {
				sw = 1;
				k  = i1;
			}
		}
	}

	return k;
}

　　提示したプログラムは，以下に示す 3 つの 2 変数関数の最小値を求めるためのものです．一様交叉，一般的な突然変異，及び，エリート＋ルーレット選択を使用しています．

  1. f = (x - 1)2 + (y - 2)2　　最小値：(1,2) で 0.0
  2. f = 100(y - x2)2 + (1 - x)2　　最小値：(1,1) で 0.0
  3. f = (1.5 - x(1 - y))2 + (2.25 - x(1 - y2))2 + (2.625 - x(1 - y3))2　　最小値：(3,0.5) で 0.0

コンパイルした後，入力データ記述ファイル「i_data」を適当に修正し，コマンドラインから，

test < i_data

と入力してやることによって実行できます．

　　入力データ記述ファイル「i_data」は，例えば，以下のように記述します． 
世代交代数 1000 乱数の初期値 12345 関数 1 出力 100
個体数 20 子供の数 20 遺伝子長 20
エリート 2 突然変異率 0.03
　　上に示した各データにおいて，日本語の部分は次に続くデータの説明になっていますので，数字の部分だけを修正してください．日本語の部分を変更しても構いませんが，削除したり，間に半角のスペースを入れるようなことはしないでください．各データの意味は以下に示す通りです．