変更履歴

240822 多峰性関数であるRastrigin function の探索結果を図5に追加しました

焼きなまし法により三次元関数の最小値を探索してみました。
使用したベンチマーク関数は粒子群最適化の場合と同じ下記の４つです。

Sphere function
Five well potential function
Himmelblau's function
Rosenbrock function

掲載した４本のプログラムは関数部分（function関数およびset_function_conditions関数）が違うだけで、他は全て同じです。
図１～図４に各関数での探索結果を示しました。
焼きなまし法は生物を模したものではないので無機的でちょっと動作がモッサリしているようにも感じます。粒子群最適化のような多点探索ではなく単点探索（ある一つの解を改良することを基本として探索をすすめる）なので、Himmelblau's function では探索結果が４か所に分かれるのが面白いところですね。多峰性関数の探索性能は、粒子群最適化よりも高いように思います。
出力結果については全部貼ると膨大な量になってしまうので、Sphere function の最初の３回分だけを貼りました。
誤りの指摘や感想等ありましたら、書き込んでいただけると幸いです。

焼きなまし法のアルゴリズムを参考にさせていただいたサイト

simulated annealing algorithm

参考にさせていただいた本

「メタヒューリスティクスとナチュラルコンピューティング」

参考にさせていただベンチマーク関数のサイト

参考にさせていただいたグラフ表示に関するサイト（実際に参考にしたのは「Visual C++の易しい使い方（23） ―三次元グラフの陰線処理―」の方ですが、今は公開されていないようです）

隠線処理を用いた三次元グラフの作成

図1　Sphere function の探索結果

図2　Five-well potential function の探索結果

図3　Himmelblau's function の探索結果

図4　Rosenbrock function の探索結果

図5 Rastrigin function の探索結果

sphere_sa.cpp

#define _USE_MATH_DEFINES
#include <math.h>
#include <iostream>
#include <fstream>
using namespace std;

const int ITERATION = 1000; // 試行回数

const int SWARM_SIZE = 20; // 群の大きさ

double X_MAX = 10; // graphのX軸の表示範囲の上限
double X_MIN = -10; // graphのX軸の表示範囲の下限

double Y_MAX = 10; // giterationaphのY軸の表示範囲の上限
double Y_MIN = -10; // graphのY軸の表示範囲の下限

double Z_MAX = 200; //  graphのZ軸の表示範囲の上限
double Z_MIN = 0; // graphのZ軸の表示範囲の下限

string str_formula = "x * x + y * y"; // グラフに表示する数式

bool b_auto_cooling_rate = true;
double INITIAL_T = 10000;
const int REPETITION = 100;
double COOLING_RATE = 0.998;
double* x_sa = nullptr;
double* y_sa = nullptr;
double* xm = nullptr;
double* ym = nullptr;
double* minim = nullptr;
double* z_sa = nullptr;
double d;
double T = 0;
int change_counter = 0;

double x[ITERATION * REPETITION][SWARM_SIZE];
double y[ITERATION * REPETITION][SWARM_SIZE];
double z[ITERATION * REPETITION][SWARM_SIZE];

double function(double x, double y)
{
	double z;

	z = x * x + y * y;

	return z;
}

void set_function_conditions()
{
	X_MAX = 10;
	X_MIN = -10;
	Y_MAX = 10;
	Y_MIN = -10;
	Z_MAX = 200;
	Z_MIN = 0;
	str_formula = "x * x + y * y";
}

double max1(double a, double b)
{
	double c = a;
	if (b > a)
		c = b;
	return c;
}

double min1(double a, double b)
{
	double c = a;
	if (b < a)
		c = b;
	return c;
}

double accept(double z, double minim, double T, double d)
{
	double p = -(z - minim) / (d * T);
	return pow(exp(1),  p);
}

double fRand(double fMin, double fMax)
{
	double f = (double)rand() / RAND_MAX;
	return fMin + f * (fMax - fMin);
}

void delete_simulated_annealing()
{
	delete []x_sa;
	delete []y_sa;
	delete []xm;
	delete []ym;
	delete []minim;
	delete []z_sa;

	x_sa = nullptr;
	y_sa = nullptr;
	xm = nullptr;
	ym = nullptr;
	minim = nullptr;
	z_sa = nullptr;
}

double get_particle_x(int i)
{
	return xm[i];
}

double get_particle_y(int i)
{
	return ym[i];
}

double get_particle_z(int i)
{
	return minim[i];
}

void init_simulated_annealing()
{
	x_sa = new double[SWARM_SIZE];
	y_sa = new double[SWARM_SIZE];
	xm = new double[SWARM_SIZE];
	ym = new double[SWARM_SIZE];
	minim = new double[SWARM_SIZE];
	z_sa = new double[SWARM_SIZE];

//	x_sa = fRand(-30, 30);
//	y_sa = fRand(-30, 30);
	for (int i = 0; i < SWARM_SIZE; i++) {
		x_sa[i] = fRand(X_MIN, X_MAX);
		y_sa[i] = fRand(Y_MIN, Y_MAX);
		xm[i] = x_sa[i];
		ym[i] = y_sa[i];
		minim[i] = function(x_sa[i], y_sa[i]);
	}
	d = (1.6 * (pow(10, -23)));
	T = INITIAL_T;

	if (b_auto_cooling_rate == true) {
		double area_scale = sqrt((X_MAX - X_MIN) * (Y_MAX - Y_MIN));
		INITIAL_T = area_scale * area_scale;
		COOLING_RATE = pow(pow((double)area_scale, 2) / (100.0 * INITIAL_T), 1.0 / ITERATION);
		T = INITIAL_T;
	}
	
	cout << "INITIAL_T = " << INITIAL_T << " COOLING_RATE = " << COOLING_RATE << " ITERATION = " << ITERATION << endl;
}

bool optimize_simulated_annealing(int i)
{
	bool changed = false;
	for (int j = 0; j < SWARM_SIZE; j++) {
//		double step = (double)X_MAX / 10;
		double x_step = (double)X_MAX * (ITERATION - i) / ITERATION;
		double y_step = (double)Y_MAX * (ITERATION - i) / ITERATION;
		x_sa[j] = x_sa[j] + fRand(-x_step, x_step);
		y_sa[j] = y_sa[j] + fRand(-y_step, y_step);
		x_sa[j] = max1(min1(x_sa[j], X_MAX), X_MIN);
		y_sa[j] = max1(min1(y_sa[j], Y_MAX), Y_MIN);
		z_sa[j] = function(x_sa[j], y_sa[j]);
		if (z_sa[j] < minim[j] || (accept(z_sa[j], minim[j], T, d) > (fRand(0, 1)))) {
//		if (z_sa[j] < minim[j]) { // Hill Climb
			changed = true;
			minim[j] = z_sa[j];
			xm[j] = x_sa[j];
			ym[j] = y_sa[j];
		}
	}
//	T = T * COOLING_RATE;

	return changed;
}

void sort_simulated_annealing()
{
	for (int i = 0; i < SWARM_SIZE - 1; i++) {
		for (int j = i + 1; j < SWARM_SIZE; j++) {
			if (minim[i] > minim[j]) {
				double tmp_x = xm[i];
				double tmp_y = ym[i];
				double tmp_z = minim[i];
				xm[i] = xm[j];
				ym[i] = ym[j];
				minim[i] = minim[j];
				xm[j] = tmp_x;
				ym[j] = tmp_y;
				minim[j] = tmp_z;
			}
		}
	}
}

void view_simulated_annealing()
{
	for (int j = 0; j < SWARM_SIZE; j++) {
		double xx = xm[j];
		double yy = ym[j];
		double zz = minim[j];
		std::cout << "x:" << xx << " y:" << yy << " z:" << zz << std::endl;
	}
}

void create_optimize_data()
{
	// 初期値の保存
	for (int i = 0; i < SWARM_SIZE; i++) {
		x[0][i] = get_particle_x(i);
		y[0][i] = get_particle_y(i);
		z[0][i] = get_particle_z(i);
	}
	change_counter = 1;

	for (int i = 0; i < ITERATION; i++) {
		for (int j = 0; j < REPETITION; j++) {
			bool changed = optimize_simulated_annealing(i);
			if (changed == true) {
				std::cout << std::endl << i << " T:" << T << std::endl;
				view_simulated_annealing();
				for (int k = 0; k < SWARM_SIZE; k++) {
					x[change_counter][k] = get_particle_x(k);
					y[change_counter][k] = get_particle_y(k);
					z[change_counter][k] = get_particle_z(k);
				}
				change_counter++;
			}
		}
		T = T * COOLING_RATE;
	}
	cout << endl;
}

void write_data()
{
	ofstream fout;
	fout.open("graph_data.txt");

	fout << str_formula << endl;
	cout << str_formula << endl ;

//	fout << ITERATION << endl;
//	cout << ITERATION << endl;
	fout << "ITERATION " << change_counter << endl;
	cout << "ITERATION " << change_counter << endl;

	fout << "SWARM_SIZE " << SWARM_SIZE << endl;
	cout << "SWARM_SIZE " << SWARM_SIZE << endl;

	fout << "X_MAX " << X_MAX << endl;
	cout << "X_MAX " << X_MAX << endl;

	fout << "X_MIN " << X_MIN << endl;
	cout << "X_MIN " << X_MIN << endl;

	fout << "Y_MAX " << Y_MAX << endl;
	cout << "Y_MAX " << Y_MAX << endl;

	fout << "Y_MIN " << Y_MIN << endl;
	cout << "Y_MIN " << Y_MIN << endl;

	fout << "Z_MAX " << Z_MAX << endl;
	cout << "Z_MAX " << Z_MAX << endl;

	fout << "Z_MIN " << Z_MIN << endl;
	cout << "Z_MIN " << Z_MIN << endl;

	fout << "DATA" << endl;
	cout << "DATA" << endl;

//	for (int i = 0; i < ITERATION; i++) {
	for (int i = 0; i < change_counter; i++) {
		cout << i << endl;
		for (int j = 0; j < SWARM_SIZE; j++) {
			cout << x[i][j] << " " << y[i][j] << " " << z[i][j] << endl;
			fout << x[i][j] << " " << y[i][j] << " " << z[i][j] << endl;
		}
		cout << endl;
		fout << endl;
	}

	fout.close();
}

void print_sorted_result()
{
	sort_simulated_annealing();
//	std::cout << std::endl;
	std::cout << "sorted result" << std::endl;
	view_simulated_annealing();
}

int main()
{
	srand ((unsigned)time(NULL));

	set_function_conditions();

	init_simulated_annealing();

	create_optimize_data();

	write_data();

	print_sorted_result();

	delete_simulated_annealing();

	return 0;
}

five_well_potential_sa.cpp

#define _USE_MATH_DEFINES
#include <math.h>
#include <iostream>
#include <fstream>
using namespace std;

const int ITERATION = 1000; // 試行回数

const int SWARM_SIZE = 20; // 群の大きさ

double X_MAX = 10; // graphのX軸の表示範囲の上限
double X_MIN = -10; // graphのX軸の表示範囲の下限

double Y_MAX = 10; // giterationaphのY軸の表示範囲の上限
double Y_MIN = -10; // graphのY軸の表示範囲の下限

double Z_MAX = 200; //  graphのZ軸の表示範囲の上限
double Z_MIN = 0; // graphのZ軸の表示範囲の下限

string str_formula = "x * x + y * y"; // グラフに表示する数式

bool b_auto_cooling_rate = true;
double INITIAL_T = 10000;
const int REPETITION = 100;
double COOLING_RATE = 0.998;
double* x_sa = nullptr;
double* y_sa = nullptr;
double* xm = nullptr;
double* ym = nullptr;
double* minim = nullptr;
double* z_sa = nullptr;
double d;
double T = 0;
int change_counter = 0;

double x[ITERATION * REPETITION][SWARM_SIZE];
double y[ITERATION * REPETITION][SWARM_SIZE];
double z[ITERATION * REPETITION][SWARM_SIZE];

double function(double x, double y)
{
	double z;

	z = (1.0 - 1.0 / (1.0 + 0.05 * (x * x + (y - 10) * (y - 10)))
		- 1.0 / (1.0 + 0.05 * ((x - 10) * (x - 10) + y * y))
		- 1.5 / (1.0 + 0.03 * ((x + 10) * (x + 10) + y * y))
		- 2.0 / (1.0 + 0.05 * ((x - 5) * (x - 5) + (y + 10) * (y + 10)))
		- 1.0 / (1.0 + 0.1 * ((x + 5) * (x + 5) + (y + 10) * (y + 10))))
		* (1.0 + 0.0001 * exp(log(x * x + y * y) * 1.2));

	return z;
}

void set_function_conditions()
{
	X_MAX = 20;
	X_MIN = -20;
	Y_MAX = 20;
	Y_MIN = -20;
	Z_MAX = 1.5;
	Z_MIN = -1.5;
	str_formula = "(1 - 1 / (1 + 0.05 * (x**2 + (y - 10)**2)) - 1 / (1 + 0.05 * ((x - 10)**2 + y**2)) - 1.5 / (1 + 0.03 * ((x + 10)**2 + y**2)) - 2 / (1 + 0.05 * ((x - 5)**2 + (y + 10)**2)) - 1 / (1 + 0.1 * ((x + 5)**2 + (y + 10)**2))) * (1 + 0.0001 * exp(log(x**2 + y**2) * 1.2))";
}

double max1(double a, double b)
{
	double c = a;
	if (b > a)
		c = b;
	return c;
}

double min1(double a, double b)
{
	double c = a;
	if (b < a)
		c = b;
	return c;
}

double accept(double z, double minim, double T, double d)
{
	double p = -(z - minim) / (d * T);
	return pow(exp(1),  p);
}

double fRand(double fMin, double fMax)
{
	double f = (double)rand() / RAND_MAX;
	return fMin + f * (fMax - fMin);
}

void delete_simulated_annealing()
{
	delete []x_sa;
	delete []y_sa;
	delete []xm;
	delete []ym;
	delete []minim;
	delete []z_sa;

	x_sa = nullptr;
	y_sa = nullptr;
	xm = nullptr;
	ym = nullptr;
	minim = nullptr;
	z_sa = nullptr;
}

double get_particle_x(int i)
{
	return xm[i];
}

double get_particle_y(int i)
{
	return ym[i];
}

double get_particle_z(int i)
{
	return minim[i];
}

void init_simulated_annealing()
{
	x_sa = new double[SWARM_SIZE];
	y_sa = new double[SWARM_SIZE];
	xm = new double[SWARM_SIZE];
	ym = new double[SWARM_SIZE];
	minim = new double[SWARM_SIZE];
	z_sa = new double[SWARM_SIZE];

//	x_sa = fRand(-30, 30);
//	y_sa = fRand(-30, 30);
	for (int i = 0; i < SWARM_SIZE; i++) {
		x_sa[i] = fRand(X_MIN, X_MAX);
		y_sa[i] = fRand(Y_MIN, Y_MAX);
		xm[i] = x_sa[i];
		ym[i] = y_sa[i];
		minim[i] = function(x_sa[i], y_sa[i]);
	}
	d = (1.6 * (pow(10, -23)));
	T = INITIAL_T;

	if (b_auto_cooling_rate == true) {
		double area_scale = sqrt((X_MAX - X_MIN) * (Y_MAX - Y_MIN));
		INITIAL_T = area_scale * area_scale;
		COOLING_RATE = pow(pow((double)area_scale, 2) / (100.0 * INITIAL_T), 1.0 / ITERATION);
		T = INITIAL_T;
	}
	
	cout << "INITIAL_T = " << INITIAL_T << " COOLING_RATE = " << COOLING_RATE << " ITERATION = " << ITERATION << endl;
}

bool optimize_simulated_annealing(int i)
{
	bool changed = false;
	for (int j = 0; j < SWARM_SIZE; j++) {
//		double step = (double)X_MAX / 10;
		double x_step = (double)X_MAX * (ITERATION - i) / ITERATION;
		double y_step = (double)Y_MAX * (ITERATION - i) / ITERATION;
		x_sa[j] = x_sa[j] + fRand(-x_step, x_step);
		y_sa[j] = y_sa[j] + fRand(-y_step, y_step);
		x_sa[j] = max1(min1(x_sa[j], X_MAX), X_MIN);
		y_sa[j] = max1(min1(y_sa[j], Y_MAX), Y_MIN);
		z_sa[j] = function(x_sa[j], y_sa[j]);
		if (z_sa[j] < minim[j] || (accept(z_sa[j], minim[j], T, d) > (fRand(0, 1)))) {
//		if (z_sa[j] < minim[j]) { // Hill Climb
			changed = true;
			minim[j] = z_sa[j];
			xm[j] = x_sa[j];
			ym[j] = y_sa[j];
		}
	}
//	T = T * COOLING_RATE;

	return changed;
}

void sort_simulated_annealing()
{
	for (int i = 0; i < SWARM_SIZE - 1; i++) {
		for (int j = i + 1; j < SWARM_SIZE; j++) {
			if (minim[i] > minim[j]) {
				double tmp_x = xm[i];
				double tmp_y = ym[i];
				double tmp_z = minim[i];
				xm[i] = xm[j];
				ym[i] = ym[j];
				minim[i] = minim[j];
				xm[j] = tmp_x;
				ym[j] = tmp_y;
				minim[j] = tmp_z;
			}
		}
	}
}

void view_simulated_annealing()
{
	for (int j = 0; j < SWARM_SIZE; j++) {
		double xx = xm[j];
		double yy = ym[j];
		double zz = minim[j];
		std::cout << "x:" << xx << " y:" << yy << " z:" << zz << std::endl;
	}
}

void create_optimize_data()
{
	// 初期値の保存
	for (int i = 0; i < SWARM_SIZE; i++) {
		x[0][i] = get_particle_x(i);
		y[0][i] = get_particle_y(i);
		z[0][i] = get_particle_z(i);
	}
	change_counter = 1;

	for (int i = 0; i < ITERATION; i++) {
		for (int j = 0; j < REPETITION; j++) {
			bool changed = optimize_simulated_annealing(i);
			if (changed == true) {
				std::cout << std::endl << i << " T:" << T << std::endl;
				view_simulated_annealing();
				for (int k = 0; k < SWARM_SIZE; k++) {
					x[change_counter][k] = get_particle_x(k);
					y[change_counter][k] = get_particle_y(k);
					z[change_counter][k] = get_particle_z(k);
				}
				change_counter++;
			}
		}
		T = T * COOLING_RATE;
	}
	cout << endl;
}

void write_data()
{
	ofstream fout;
	fout.open("graph_data.txt");

	fout << str_formula << endl;
	cout << str_formula << endl ;

//	fout << ITERATION << endl;
//	cout << ITERATION << endl;
	fout << "ITERATION " << change_counter << endl;
	cout << "ITERATION " << change_counter << endl;

	fout << "SWARM_SIZE " << SWARM_SIZE << endl;
	cout << "SWARM_SIZE " << SWARM_SIZE << endl;

	fout << "X_MAX " << X_MAX << endl;
	cout << "X_MAX " << X_MAX << endl;

	fout << "X_MIN " << X_MIN << endl;
	cout << "X_MIN " << X_MIN << endl;

	fout << "Y_MAX " << Y_MAX << endl;
	cout << "Y_MAX " << Y_MAX << endl;

	fout << "Y_MIN " << Y_MIN << endl;
	cout << "Y_MIN " << Y_MIN << endl;

	fout << "Z_MAX " << Z_MAX << endl;
	cout << "Z_MAX " << Z_MAX << endl;

	fout << "Z_MIN " << Z_MIN << endl;
	cout << "Z_MIN " << Z_MIN << endl;

	fout << "DATA" << endl;
	cout << "DATA" << endl;

//	for (int i = 0; i < ITERATION; i++) {
	for (int i = 0; i < change_counter; i++) {
		cout << i << endl;
		for (int j = 0; j < SWARM_SIZE; j++) {
			cout << x[i][j] << " " << y[i][j] << " " << z[i][j] << endl;
			fout << x[i][j] << " " << y[i][j] << " " << z[i][j] << endl;
		}
		cout << endl;
		fout << endl;
	}

	fout.close();
}

void print_sorted_result()
{
	sort_simulated_annealing();
//	std::cout << std::endl;
	std::cout << "sorted result" << std::endl;
	view_simulated_annealing();
}

int main()
{
	srand ((unsigned)time(NULL));

	set_function_conditions();

	init_simulated_annealing();

	create_optimize_data();

	write_data();

	print_sorted_result();

	delete_simulated_annealing();

	return 0;
}

himmelblau_sa.cpp

#define _USE_MATH_DEFINES
#include <math.h>
#include <iostream>
#include <fstream>
using namespace std;

const int ITERATION = 1000; // 試行回数

const int SWARM_SIZE = 20; // 群の大きさ

double X_MAX = 10; // graphのX軸の表示範囲の上限
double X_MIN = -10; // graphのX軸の表示範囲の下限

double Y_MAX = 10; // giterationaphのY軸の表示範囲の上限
double Y_MIN = -10; // graphのY軸の表示範囲の下限

double Z_MAX = 200; //  graphのZ軸の表示範囲の上限
double Z_MIN = 0; // graphのZ軸の表示範囲の下限

string str_formula = "x * x + y * y"; // グラフに表示する数式

bool b_auto_cooling_rate = true;
double INITIAL_T = 10000;
const int REPETITION = 100;
double COOLING_RATE = 0.998;
double* x_sa = nullptr;
double* y_sa = nullptr;
double* xm = nullptr;
double* ym = nullptr;
double* minim = nullptr;
double* z_sa = nullptr;
double d;
double T = 0;
int change_counter = 0;

double x[ITERATION * REPETITION][SWARM_SIZE];
double y[ITERATION * REPETITION][SWARM_SIZE];
double z[ITERATION * REPETITION][SWARM_SIZE];

double function(double x, double y)
{
	double z;

	z = pow((x * x + y - 11), 2) + pow((x + y * y - 7), 2);

	return z;
}

void set_function_conditions()
{
	X_MAX = 5;
	X_MIN = -5;
	Y_MAX = 5;
	Y_MIN = -5;
	Z_MAX = 500;
	Z_MIN = 0;
	str_formula = "(x**2 + y - 11)**2 + (x + y**2 - 7)**2";
}

double max1(double a, double b)
{
	double c = a;
	if (b > a)
		c = b;
	return c;
}

double min1(double a, double b)
{
	double c = a;
	if (b < a)
		c = b;
	return c;
}

double accept(double z, double minim, double T, double d)
{
	double p = -(z - minim) / (d * T);
	return pow(exp(1),  p);
}

double fRand(double fMin, double fMax)
{
	double f = (double)rand() / RAND_MAX;
	return fMin + f * (fMax - fMin);
}

void delete_simulated_annealing()
{
	delete []x_sa;
	delete []y_sa;
	delete []xm;
	delete []ym;
	delete []minim;
	delete []z_sa;

	x_sa = nullptr;
	y_sa = nullptr;
	xm = nullptr;
	ym = nullptr;
	minim = nullptr;
	z_sa = nullptr;
}

double get_particle_x(int i)
{
	return xm[i];
}

double get_particle_y(int i)
{
	return ym[i];
}

double get_particle_z(int i)
{
	return minim[i];
}

void init_simulated_annealing()
{
	x_sa = new double[SWARM_SIZE];
	y_sa = new double[SWARM_SIZE];
	xm = new double[SWARM_SIZE];
	ym = new double[SWARM_SIZE];
	minim = new double[SWARM_SIZE];
	z_sa = new double[SWARM_SIZE];

//	x_sa = fRand(-30, 30);
//	y_sa = fRand(-30, 30);
	for (int i = 0; i < SWARM_SIZE; i++) {
		x_sa[i] = fRand(X_MIN, X_MAX);
		y_sa[i] = fRand(Y_MIN, Y_MAX);
		xm[i] = x_sa[i];
		ym[i] = y_sa[i];
		minim[i] = function(x_sa[i], y_sa[i]);
	}
	d = (1.6 * (pow(10, -23)));
	T = INITIAL_T;

	if (b_auto_cooling_rate == true) {
		double area_scale = sqrt((X_MAX - X_MIN) * (Y_MAX - Y_MIN));
		INITIAL_T = area_scale * area_scale;
		COOLING_RATE = pow(pow((double)area_scale, 2) / (100.0 * INITIAL_T), 1.0 / ITERATION);
		T = INITIAL_T;
	}
	
	cout << "INITIAL_T = " << INITIAL_T << " COOLING_RATE = " << COOLING_RATE << " ITERATION = " << ITERATION << endl;
}

bool optimize_simulated_annealing(int i)
{
	bool changed = false;
	for (int j = 0; j < SWARM_SIZE; j++) {
//		double step = (double)X_MAX / 10;
		double x_step = (double)X_MAX * (ITERATION - i) / ITERATION;
		double y_step = (double)Y_MAX * (ITERATION - i) / ITERATION;
		x_sa[j] = x_sa[j] + fRand(-x_step, x_step);
		y_sa[j] = y_sa[j] + fRand(-y_step, y_step);
		x_sa[j] = max1(min1(x_sa[j], X_MAX), X_MIN);
		y_sa[j] = max1(min1(y_sa[j], Y_MAX), Y_MIN);
		z_sa[j] = function(x_sa[j], y_sa[j]);
		if (z_sa[j] < minim[j] || (accept(z_sa[j], minim[j], T, d) > (fRand(0, 1)))) {
//		if (z_sa[j] < minim[j]) { // Hill Climb
			changed = true;
			minim[j] = z_sa[j];
			xm[j] = x_sa[j];
			ym[j] = y_sa[j];
		}
	}
//	T = T * COOLING_RATE;

	return changed;
}

void sort_simulated_annealing()
{
	for (int i = 0; i < SWARM_SIZE - 1; i++) {
		for (int j = i + 1; j < SWARM_SIZE; j++) {
			if (minim[i] > minim[j]) {
				double tmp_x = xm[i];
				double tmp_y = ym[i];
				double tmp_z = minim[i];
				xm[i] = xm[j];
				ym[i] = ym[j];
				minim[i] = minim[j];
				xm[j] = tmp_x;
				ym[j] = tmp_y;
				minim[j] = tmp_z;
			}
		}
	}
}

void view_simulated_annealing()
{
	for (int j = 0; j < SWARM_SIZE; j++) {
		double xx = xm[j];
		double yy = ym[j];
		double zz = minim[j];
		std::cout << "x:" << xx << " y:" << yy << " z:" << zz << std::endl;
	}
}

void create_optimize_data()
{
	// 初期値の保存
	for (int i = 0; i < SWARM_SIZE; i++) {
		x[0][i] = get_particle_x(i);
		y[0][i] = get_particle_y(i);
		z[0][i] = get_particle_z(i);
	}
	change_counter = 1;

	for (int i = 0; i < ITERATION; i++) {
		for (int j = 0; j < REPETITION; j++) {
			bool changed = optimize_simulated_annealing(i);
			if (changed == true) {
				std::cout << std::endl << i << " T:" << T << std::endl;
				view_simulated_annealing();
				for (int k = 0; k < SWARM_SIZE; k++) {
					x[change_counter][k] = get_particle_x(k);
					y[change_counter][k] = get_particle_y(k);
					z[change_counter][k] = get_particle_z(k);
				}
				change_counter++;
			}
		}
		T = T * COOLING_RATE;
	}
	cout << endl;
}

void write_data()
{
	ofstream fout;
	fout.open("graph_data.txt");

	fout << str_formula << endl;
	cout << str_formula << endl ;

//	fout << ITERATION << endl;
//	cout << ITERATION << endl;
	fout << "ITERATION " << change_counter << endl;
	cout << "ITERATION " << change_counter << endl;

	fout << "SWARM_SIZE " << SWARM_SIZE << endl;
	cout << "SWARM_SIZE " << SWARM_SIZE << endl;

	fout << "X_MAX " << X_MAX << endl;
	cout << "X_MAX " << X_MAX << endl;

	fout << "X_MIN " << X_MIN << endl;
	cout << "X_MIN " << X_MIN << endl;

	fout << "Y_MAX " << Y_MAX << endl;
	cout << "Y_MAX " << Y_MAX << endl;

	fout << "Y_MIN " << Y_MIN << endl;
	cout << "Y_MIN " << Y_MIN << endl;

	fout << "Z_MAX " << Z_MAX << endl;
	cout << "Z_MAX " << Z_MAX << endl;

	fout << "Z_MIN " << Z_MIN << endl;
	cout << "Z_MIN " << Z_MIN << endl;

	fout << "DATA" << endl;
	cout << "DATA" << endl;

//	for (int i = 0; i < ITERATION; i++) {
	for (int i = 0; i < change_counter; i++) {
		cout << i << endl;
		for (int j = 0; j < SWARM_SIZE; j++) {
			cout << x[i][j] << " " << y[i][j] << " " << z[i][j] << endl;
			fout << x[i][j] << " " << y[i][j] << " " << z[i][j] << endl;
		}
		cout << endl;
		fout << endl;
	}

	fout.close();
}

void print_sorted_result()
{
	sort_simulated_annealing();
//	std::cout << std::endl;
	std::cout << "sorted result" << std::endl;
	view_simulated_annealing();
}

int main()
{
	srand ((unsigned)time(NULL));

	set_function_conditions();

	init_simulated_annealing();

	create_optimize_data();

	write_data();

	print_sorted_result();

	delete_simulated_annealing();

	return 0;
}

rosenbrock_sa.cpp

#define _USE_MATH_DEFINES
#include <math.h>
#include <iostream>
#include <fstream>
using namespace std;

const int ITERATION = 1000; // 試行回数

const int SWARM_SIZE = 20; // 群の大きさ

double X_MAX = 10; // graphのX軸の表示範囲の上限
double X_MIN = -10; // graphのX軸の表示範囲の下限

double Y_MAX = 10; // giterationaphのY軸の表示範囲の上限
double Y_MIN = -10; // graphのY軸の表示範囲の下限

double Z_MAX = 200; //  graphのZ軸の表示範囲の上限
double Z_MIN = 0; // graphのZ軸の表示範囲の下限

string str_formula = "x * x + y * y"; // グラフに表示する数式

double x[ITERATION][SWARM_SIZE];
double y[ITERATION][SWARM_SIZE];
double z[ITERATION][SWARM_SIZE];

double function(double x, double y)
{
	double z;

	z = 100 * (y - x * x) * (y  - x * x) + (1 - x) * (1 - x);

	return z;
}

void set_function_conditions()
{
	X_MAX = 2.048;
	X_MIN = -2.048;
	Y_MAX = 2.048;
	Y_MIN = -2.048;
	Z_MAX = 4000;
	Z_MIN = 0;
	str_formula = "100 * (y - x * x)**2 + (1 - x)**2";
}

int N_P = SWARM_SIZE * 2; // number of population

int CUT_LINE = SWARM_SIZE;
double A_GA = 0.3;
double B_GA = 0.1;
double MUTATION_ODDS = 0.1;
int change_counter = 0;

double get_random()
{
	return (double)rand() / ((double)RAND_MAX + 1);
}

double max1(double a, double b)
{
	double c = a;
	if (b > a)
		c = b;
	return c;
}

double min1(double a, double b)
{
	double c = a;
	if (b < a)
		c = b;
	return c;
}

class POPULATION
{
protected:
	double x, y, z;
	int* gene;
	long long weight_p;
	long long value_p;
public:
	POPULATION();
	~POPULATION();
	void calc_weight_and_value();
	long long get_value() { return value_p; }
	long long get_weight() { return weight_p; }
	void set_gene(int i, int g);
	int get_gene(int i) { return gene[i]; }
	void mutation();
	void initialize();
	void set_xyz(double x1, double y1, double z1) { x = x1; y = y1; z = z1; }
	double get_x() { return x; }
	double get_y() { return y; }
	double get_z() { return z; }
};

POPULATION** population = nullptr;

POPULATION::POPULATION()
{

}

POPULATION::~POPULATION()
{

}

void POPULATION::initialize()
{
	double x1 = get_random() * (X_MAX - X_MIN) + X_MIN;
	double y1 = get_random() * (Y_MAX - Y_MIN) + Y_MIN;
	double z1 = function(x1, y1);
	set_xyz(x1, y1, z1);
}

void POPULATION::mutation()
{
	initialize();
}

void create_object()
{
	population = new POPULATION*[N_P];
	for (int i = 0; i < N_P; i++)
		population[i] = new POPULATION;
}

void initialize_population()
{
	for (int i = 0; i < N_P; i++) 
		population[i]->initialize();
}

void ranking()
{
	for (int i = 0; i < N_P - 1; i++) {
		for (int j = i + 1; j < N_P; j++) {
			if (population[i]->get_z() > population[j]->get_z()) {
				POPULATION* tmp = population[i];
				population[i] = population[j];
				population[j] = tmp;
			}
		}
	}
}

void cross_over(int i)
{
	for (int j = CUT_LINE; j < N_P; j++) {
		int parent1 = rand() % CUT_LINE;
		int parent2 = rand() % CUT_LINE;
	
		double x1 = population[parent1]->get_x();
		double y1 = population[parent1]->get_y();
		double x2 = population[parent2]->get_x();
		double y2 = population[parent2]->get_y();
		double dx = abs(x1 - x2);
		double dy = abs(y1 - y2);
		double min_x = min(x1, x2);
		double min_y = min(y1, y2);
		double max_x = max(x1, x2);
		double max_y = max(y1, y2);
		double min_cx = min_x - A_GA * dx;
		double max_cx = max_x + A_GA * dx;
		double min_cy = min_y - A_GA * dy;
		double max_cy = max_y + A_GA * dy;
//		double next_x = get_random() * (max_cx - min_cx) + min_cx;
//		double next_y = get_random() * (max_cy - min_cy) + min_cy;
		double DX = (X_MAX - X_MIN) * B_GA * (get_random() - 0.5);
		double DY = (Y_MAX - Y_MIN) * B_GA * (get_random() - 0.5);
		double next_x = get_random() * (max_cx - min_cx) + min_cx + DX;
		double next_y = get_random() * (max_cy - min_cy) + min_cy + DY;
		next_x = max1(min1(next_x, X_MAX), X_MIN);
		next_y = max1(min1(next_y, Y_MAX), Y_MIN);
		double next_z = function(next_x, next_y);
		population[j]->set_xyz(next_x, next_y, next_z);

		if (get_random() < MUTATION_ODDS)
			population[j]->mutation();
	}
}

void init_genetic_algorithm()
{
	create_object();
	initialize_population();
	ranking();
	change_counter = 0;
	for (int j = 0; j < CUT_LINE; j++) {
		x[change_counter][j] = population[j]->get_x();
		y[change_counter][j] = population[j]->get_y();
		z[change_counter][j] = population[j]->get_z();
	}
	change_counter++;		
}

void delete_genetic_algorithm()
{
	for (int i = 0; i < N_P; i++)
		delete population[i];
	delete []population;
	population = nullptr;
}

void save_result(int i)
{
//	std::cout << std::endl << i << std::endl;
	bool b_changed = false;
	for (int j = 0; j < CUT_LINE; j++) {
		double xx = population[j]->get_x();
		double yy = population[j]->get_y();
		double zz = population[j]->get_z();
//		std::cout << "x:" << xx << " y:" << yy << " z:" << zz << std::endl;

		if ((float)xx != (float)x[change_counter - 1][j])
			b_changed = true;
 		if ((float)yy != (float)y[change_counter - 1][j])
			b_changed = true;
		if ((float)zz != (float)z[change_counter - 1][j])
			b_changed = true;
	}

	if (b_changed == true) {
		for (int j = 0; j < CUT_LINE; j++) {
			x[change_counter][j] = population[j]->get_x();
			y[change_counter][j] = population[j]->get_y();
			z[change_counter][j] = population[j]->get_z();
		}
		change_counter++;
	}		
}

void optimize_genetic_algorithm()
{
	for (int i = 0; i < ITERATION; i++) {
		cross_over(i);
		ranking();
		save_result(i);
	}
}	

void write_data();

int main()
{
	srand ((unsigned)time(NULL));

	set_function_conditions();

	init_genetic_algorithm();

	optimize_genetic_algorithm();

	write_data();

	delete_genetic_algorithm();

	return 0;
}

void write_data()
{
	ofstream fout;
	fout.open("graph_data.txt");

	fout << str_formula << endl;
	cout << str_formula << endl ;

//	fout << ITERATION << endl;
//	cout << ITERATION << endl;
	fout << "ITERATION " << change_counter << endl;
	cout << "ITERATION " << change_counter << endl;

	fout << "SWARM_SIZE " << SWARM_SIZE << endl;
	cout << "SWARM_SIZE " << SWARM_SIZE << endl;

	fout << "X_MAX " << X_MAX << endl;
	cout << "X_MAX " << X_MAX << endl;

	fout << "X_MIN " << X_MIN << endl;
	cout << "X_MIN " << X_MIN << endl;

	fout << "Y_MAX " << Y_MAX << endl;
	cout << "Y_MAX " << Y_MAX << endl;

	fout << "Y_MIN " << Y_MIN << endl;
	cout << "Y_MIN " << Y_MIN << endl;

	fout << "Z_MAX " << Z_MAX << endl;
	cout << "Z_MAX " << Z_MAX << endl;

	fout << "Z_MIN " << Z_MIN << endl;
	cout << "Z_MIN " << Z_MIN << endl;

	fout << "DATA" << endl;
	cout << "DATA" << endl;

//	for (int i = 0; i < ITERATION; i++) {
	for (int i = 0; i < change_counter; i++) {
		cout << i << endl;
		for (int j = 0; j < SWARM_SIZE; j++) {
			cout << x[i][j] << " " << y[i][j] << " " << z[i][j] << endl;
			fout << x[i][j] << " " << y[i][j] << " " << z[i][j] << endl;
		}
		cout << endl;
		fout << endl;
	}

	fout.close();
}

rastrigin_sa.cpp

#define _USE_MATH_DEFINES
#include <math.h>
#include <iostream>
#include <fstream>
using namespace std;

const int ITERATION = 1000; // 試行回数

const int SWARM_SIZE = 20; // 群の大きさ

double X_MAX = 10; // graphのX軸の表示範囲の上限
double X_MIN = -10; // graphのX軸の表示範囲の下限

double Y_MAX = 10; // giterationaphのY軸の表示範囲の上限
double Y_MIN = -10; // graphのY軸の表示範囲の下限

double Z_MAX = 200; //  graphのZ軸の表示範囲の上限
double Z_MIN = 0; // graphのZ軸の表示範囲の下限

string str_formula = "x * x + y * y"; // グラフに表示する数式

bool b_auto_cooling_rate = true;
double INITIAL_T = 10000;
const int REPETITION = 100;
double COOLING_RATE = 0.998;
double* x_sa = nullptr;
double* y_sa = nullptr;
double* xm = nullptr;
double* ym = nullptr;
double* minim = nullptr;
double* z_sa = nullptr;
double d;
double T = 0;
int change_counter = 0;

double x[ITERATION * REPETITION][SWARM_SIZE];
double y[ITERATION * REPETITION][SWARM_SIZE];
double z[ITERATION * REPETITION][SWARM_SIZE];

double function(double x, double y)
{
	double z;

	z = 20.0 + x * x - 10.0 * cos(2.0 * M_PI * x) + y * y - 10.0 * cos(2.0 * M_PI * y);

	return z;
}

void set_function_conditions()
{
	X_MAX = 10;
	X_MIN = -10;
	Y_MAX = 10;
	Y_MIN = -10;
	Z_MAX = 200;
	Z_MIN = 0;
	str_formula = "20.0 + x * x - 10.0 * cos(2.0 * pi * x) + y * y - 10.0 * cos(2.0 * pi * y)";
}

double max1(double a, double b)
{
	double c = a;
	if (b > a)
		c = b;
	return c;
}

double min1(double a, double b)
{
	double c = a;
	if (b < a)
		c = b;
	return c;
}

double accept(double z, double minim, double T, double d)
{
	double p = -(z - minim) / (d * T);
	return pow(exp(1),  p);
}

double fRand(double fMin, double fMax)
{
	double f = (double)rand() / RAND_MAX;
	return fMin + f * (fMax - fMin);
}

void delete_simulated_annealing()
{
	delete []x_sa;
	delete []y_sa;
	delete []xm;
	delete []ym;
	delete []minim;
	delete []z_sa;

	x_sa = nullptr;
	y_sa = nullptr;
	xm = nullptr;
	ym = nullptr;
	minim = nullptr;
	z_sa = nullptr;
}

double get_particle_x(int i)
{
	return xm[i];
}

double get_particle_y(int i)
{
	return ym[i];
}

double get_particle_z(int i)
{
	return minim[i];
}

void init_simulated_annealing()
{
	x_sa = new double[SWARM_SIZE];
	y_sa = new double[SWARM_SIZE];
	xm = new double[SWARM_SIZE];
	ym = new double[SWARM_SIZE];
	minim = new double[SWARM_SIZE];
	z_sa = new double[SWARM_SIZE];

//	x_sa = fRand(-30, 30);
//	y_sa = fRand(-30, 30);
	for (int i = 0; i < SWARM_SIZE; i++) {
		x_sa[i] = fRand(X_MIN, X_MAX);
		y_sa[i] = fRand(Y_MIN, Y_MAX);
		xm[i] = x_sa[i];
		ym[i] = y_sa[i];
		minim[i] = function(x_sa[i], y_sa[i]);
	}
	d = (1.6 * (pow(10, -23)));
	T = INITIAL_T;

	if (b_auto_cooling_rate == true) {
		double area_scale = sqrt((X_MAX - X_MIN) * (Y_MAX - Y_MIN));
		INITIAL_T = area_scale * area_scale;
		COOLING_RATE = pow(pow((double)area_scale, 2) / (100.0 * INITIAL_T), 1.0 / ITERATION);
		T = INITIAL_T;
	}
	
	cout << "INITIAL_T = " << INITIAL_T << " COOLING_RATE = " << COOLING_RATE << " ITERATION = " << ITERATION << endl;
}

bool optimize_simulated_annealing(int i)
{
	bool changed = false;
	for (int j = 0; j < SWARM_SIZE; j++) {
//		double step = (double)X_MAX / 10;
		double x_step = (double)X_MAX * (ITERATION - i) / ITERATION;
		double y_step = (double)Y_MAX * (ITERATION - i) / ITERATION;
		x_sa[j] = x_sa[j] + fRand(-x_step, x_step);
		y_sa[j] = y_sa[j] + fRand(-y_step, y_step);
		x_sa[j] = max1(min1(x_sa[j], X_MAX), X_MIN);
		y_sa[j] = max1(min1(y_sa[j], Y_MAX), Y_MIN);
		z_sa[j] = function(x_sa[j], y_sa[j]);
		if (z_sa[j] < minim[j] || (accept(z_sa[j], minim[j], T, d) > (fRand(0, 1)))) {
//		if (z_sa[j] < minim[j]) { // Hill Climb
			changed = true;
			minim[j] = z_sa[j];
			xm[j] = x_sa[j];
			ym[j] = y_sa[j];
		}
	}
//	T = T * COOLING_RATE;

	return changed;
}

void sort_simulated_annealing()
{
	for (int i = 0; i < SWARM_SIZE - 1; i++) {
		for (int j = i + 1; j < SWARM_SIZE; j++) {
			if (minim[i] > minim[j]) {
				double tmp_x = xm[i];
				double tmp_y = ym[i];
				double tmp_z = minim[i];
				xm[i] = xm[j];
				ym[i] = ym[j];
				minim[i] = minim[j];
				xm[j] = tmp_x;
				ym[j] = tmp_y;
				minim[j] = tmp_z;
			}
		}
	}
}

void view_simulated_annealing()
{
	for (int j = 0; j < SWARM_SIZE; j++) {
		double xx = xm[j];
		double yy = ym[j];
		double zz = minim[j];
		std::cout << "x:" << xx << " y:" << yy << " z:" << zz << std::endl;
	}
}

void create_optimize_data()
{
	// 初期値の保存
	for (int i = 0; i < SWARM_SIZE; i++) {
		x[0][i] = get_particle_x(i);
		y[0][i] = get_particle_y(i);
		z[0][i] = get_particle_z(i);
	}
	change_counter = 1;

	for (int i = 0; i < ITERATION; i++) {
		for (int j = 0; j < REPETITION; j++) {
			bool changed = optimize_simulated_annealing(i);
			if (changed == true) {
				std::cout << std::endl << i << " T:" << T << std::endl;
				view_simulated_annealing();
				for (int k = 0; k < SWARM_SIZE; k++) {
					x[change_counter][k] = get_particle_x(k);
					y[change_counter][k] = get_particle_y(k);
					z[change_counter][k] = get_particle_z(k);
				}
				change_counter++;
			}
		}
		T = T * COOLING_RATE;
	}
	cout << endl;
}

void write_data()
{
	ofstream fout;
	fout.open("graph_data.txt");

	fout << str_formula << endl;
	cout << str_formula << endl ;

//	fout << ITERATION << endl;
//	cout << ITERATION << endl;
	fout << "ITERATION " << change_counter << endl;
	cout << "ITERATION " << change_counter << endl;

	fout << "SWARM_SIZE " << SWARM_SIZE << endl;
	cout << "SWARM_SIZE " << SWARM_SIZE << endl;

	fout << "X_MAX " << X_MAX << endl;
	cout << "X_MAX " << X_MAX << endl;

	fout << "X_MIN " << X_MIN << endl;
	cout << "X_MIN " << X_MIN << endl;

	fout << "Y_MAX " << Y_MAX << endl;
	cout << "Y_MAX " << Y_MAX << endl;

	fout << "Y_MIN " << Y_MIN << endl;
	cout << "Y_MIN " << Y_MIN << endl;

	fout << "Z_MAX " << Z_MAX << endl;
	cout << "Z_MAX " << Z_MAX << endl;

	fout << "Z_MIN " << Z_MIN << endl;
	cout << "Z_MIN " << Z_MIN << endl;

	fout << "DATA" << endl;
	cout << "DATA" << endl;

//	for (int i = 0; i < ITERATION; i++) {
	for (int i = 0; i < change_counter; i++) {
		cout << i << endl;
		for (int j = 0; j < SWARM_SIZE; j++) {
			cout << x[i][j] << " " << y[i][j] << " " << z[i][j] << endl;
			fout << x[i][j] << " " << y[i][j] << " " << z[i][j] << endl;
		}
		cout << endl;
		fout << endl;
	}

	fout.close();
}

void print_sorted_result()
{
	sort_simulated_annealing();
//	std::cout << std::endl;
	std::cout << "sorted result" << std::endl;
	view_simulated_annealing();
}

int main()
{
	srand ((unsigned)time(NULL));

	set_function_conditions();

	init_simulated_annealing();

	create_optimize_data();

	write_data();

	print_sorted_result();

	delete_simulated_annealing();

	return 0;
}

出力結果 (sphere_sa)

x * x + y * y
ITERATION 177
SWARM_SIZE 20
X_MAX 10
X_MIN -10
Y_MAX 10
Y_MIN -10
Z_MAX 200
Z_MIN 0
DATA
0.323191 6.47572 42.0394
5.54064 -7.43278 85.9449
-1.53722 -9.06735 84.58
-0.199896 -5.87756 34.5857
-0.884121 -1.94494 4.56448
-3.92865 -2.25806 20.5331
-5.10178 -1.05014 27.1309
4.16974 -5.42039 46.7674
0.558794 -5.26841 28.0684
-5.33616 1.1771 29.8602
-8.27815 -7.20389 120.424
-7.43278 -9.09482 137.962
-6.95608 -5.15183 74.9285
-4.25764 -8.76034 94.871
8.90866 -2.68532 86.5752
-3.83343 -7.53166 71.4211
-0.212104 -5.23301 27.4294
-7.83013 3.4788 73.4131
-2.60903 6.97134 55.4066
9.77599 0.234077 95.6249

0.323191 6.47572 42.0394
5.54064 -7.43278 85.9449
-2.43233 -4.10352 22.7551
-0.199896 -5.87756 34.5857
-0.884121 -1.94494 4.56448
-3.92865 -2.25806 20.5331
-5.10178 -1.05014 27.1309
4.16974 -5.42039 46.7674
0.558794 -5.26841 28.0684
-4.40321 -2.98166 28.2785
0.134281 -9.80743 96.2037
-2.83456 -1.3123 9.75685
-6.95608 -5.15183 74.9285
-4.25764 -8.76034 94.871
8.90866 -2.68532 86.5752
-3.83343 -7.53166 71.4211
-0.212104 -5.23301 27.4294
-3.68175 4.07544 30.1645
-2.60903 6.97134 55.4066
9.77599 0.234077 95.6249

0.323191 6.47572 42.0394
5.54064 -7.43278 85.9449
-2.43233 -4.10352 22.7551
-0.199896 -5.87756 34.5857
-0.884121 -1.94494 4.56448
-3.53221 0.596026 12.8318
-5.10178 -1.05014 27.1309
4.16974 -5.42039 46.7674
1.46458 0.481887 2.37722
-4.40321 -2.98166 28.2785
0.134281 -9.80743 96.2037
-2.83456 -1.3123 9.75685
-6.95608 -5.15183 74.9285
1.99805 -6.54012 46.7653
8.90866 -2.68532 86.5752
-3.83343 -7.53166 71.4211
-0.212104 -5.23301 27.4294
3.19193 1.65197 12.9174
-2.60903 6.97134 55.4066
9.77599 0.234077 95.6249