変更履歴

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

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

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

掲載した４本のプログラムは関数部分（function関数およびset_function_conditions関数）が違うだけで、他は全て同じです。
図１～図４に各関数での探索結果を示しました。
粒子群最適化は単峰性の関数は得意だけれども多峰性の関数は苦手らしいですが、実際そのように感じました。（今回はあまり多峰性の関数は使用しませんでしたが）
できればリアルタイムの追従性能などについても調べてみたいです。
出力結果については全部貼ると膨大な量になってしまうので、Sphere function の最初の３回分だけを貼りました。

粒子群最適化については、下記の本を参考にさせていただきました。

「進化計算アルゴリズム入門生物の行動科学から導く最適解」

ベンチマーク関数については、下記のサイトを参考にさせていただきました。

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

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

図1 Sphere function の探索結果

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

図3 Himmelblau's function の探索結果

図4 Rosenbrock function の探索結果

図5 Rastrigin function の探索結果

sphere_ps.cpp

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

// 0以上1以下の実数乱数
#define RAND_01 ((double)rand() / RAND_MAX)

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

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 INERTIA = 0.75; // 慣性係数
double ACCEL_P = 1.0; // 加速係数（パーソナルベスト）
double ACCEL_G = 1.0; // 加速係数（グローバルベスト）

double** pos = nullptr; // 位置
double** velocity = nullptr; // 速度
double* value = nullptr; // 評価値
double** pBestPos = nullptr; // パーソナルベストの位置
double* pBestValue = nullptr; // パーソナルベストの評価値

double gBestPos[2]; // グローバルベストの位置
double gBestValue; // グローバルベストの評価値

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

// 評価値を算出し，パーソナルベストを更新する
void evaluate(int target)
{
	// 評価値を算出する
	value[target] = function(pos[target][0], pos[target][1]);

	// パーソナルベストを更新する
	if (pBestValue[target] > value[target]) {
		for (int i = 0; i < 2; i++)
			pBestPos[target][i] = pos[target][i];
		pBestValue[target] = value[target];
	}
}

void create_particle()
{
	pos = new double*[SWARM_SIZE];
	for (int i = 0; i < SWARM_SIZE; i++) {
		pos[i] = new double[2];
		pos[i][0] = X_MIN + (X_MAX - X_MIN) * RAND_01;
		pos[i][1] = Y_MIN + (Y_MAX - Y_MIN) * RAND_01;
	}

	velocity = new double*[SWARM_SIZE];
	for (int i = 0; i < SWARM_SIZE; i++) {
		velocity[i] = new double[2];
		velocity[i][0] = X_MIN + (X_MAX - X_MIN) * RAND_01;
		velocity[i][1] = Y_MIN + (Y_MAX - Y_MIN) * RAND_01;
	}

	value = new double[SWARM_SIZE];

	pBestPos = new double*[SWARM_SIZE];
	for (int i = 0; i < SWARM_SIZE; i++)
		pBestPos[i] = new double[2];

	pBestValue = new double[SWARM_SIZE];
	for (int i = 0; i < SWARM_SIZE; i++)
		pBestValue[i] = DBL_MAX;

	for (int i = 0; i < SWARM_SIZE; i++)
		evaluate(i);
}

// 個々の粒子を移動する
void move_particle(int target)
{
	for (int j = 0; j < 2; j++) {
		velocity[target][j] = INERTIA * velocity[target][j] 
			+ ACCEL_G * (gBestPos[j] - pos[target][j]) * RAND_01
			+ ACCEL_P * (pBestPos[target][j] - pos[target][j]) * RAND_01;
		pos[target][j] += velocity[target][j];
	}
	pos[target][0] = max1(min1(pos[target][0], X_MAX), X_MIN);
	pos[target][1] = max1(min1(pos[target][1], Y_MAX), Y_MIN);
	evaluate(target);
}

// 粒子群の生成
void create_particles()
{
	create_particle();

	int best = 0;
	for (int i = 0; i < SWARM_SIZE; i++)
		if (value[best] > value[i])
			best = i;

	for (int i = 0; i < 2; i++)
		gBestPos[i] = pos[best][i];

	gBestValue = value[best];
}

// 粒子群を初期化するため、粒子群を削除する
void delete_particles()
{
	for (int i = 0; i < SWARM_SIZE; i++)
		delete pos[i];
	delete []pos;
	pos = nullptr;

	for (int i = 0; i < SWARM_SIZE; i++)
		delete velocity[i];
	delete []velocity;
	velocity = nullptr;

	delete []value;
	value = nullptr;

	for (int i = 0; i < SWARM_SIZE; i++)
		delete pBestPos[i];
	delete []pBestPos;
	pBestPos = nullptr;

	delete []pBestValue;
	pBestValue = nullptr;
}

// 粒子を移動し，グローバルベストを更新する
void move_particles()
{
	// すべての粒子を移動する
	int best = -1;
	for (int j = 0; j < SWARM_SIZE; j++) {
		move_particle(j);
		if (gBestValue > value[j])
			best = j;
	}

	// グローバルベストを更新する
	if (best != -1) {
		gBestPos[0] = pos[best][0];
		gBestPos[1] = pos[best][1];
		gBestValue = value[best];
	}
}

// パーティクルを評価値の高い順にソートする
void sort_particles()
{
	for (int i = 0; i < SWARM_SIZE - 1; i++) {
		for (int j = i + 1; j < SWARM_SIZE; j++) {
			if (value[i] > value[j]) {
				double tmp_pos[2];
				double tmp_velocity[2];
				double tmp_value;
				double tmp_pBestPos[2];
				double tmp_pBestValue;

				tmp_pos[0] = pos[i][0];
				tmp_pos[1] = pos[i][1];
				tmp_velocity[0] = velocity[i][0];
				tmp_velocity[1] = velocity[i][1];
				tmp_value = value[i];
				tmp_pBestPos[0] = pBestPos[i][0];
				tmp_pBestPos[1] = pBestPos[i][1];
				tmp_pBestValue = pBestValue[i];

				pos[i][0] = pos[j][0];
				pos[i][1] = pos[j][1];
				velocity[i][0] = velocity[j][0];
				velocity[i][1] = velocity[j][1];
				value[i] = value[j];
				pBestPos[i][0] = pBestPos[j][0];
				pBestPos[i][1] = pBestPos[j][1];
				pBestValue[i] = pBestValue[j];

				pos[j][0] = tmp_pos[0];
				pos[j][1] = tmp_pos[1];
				velocity[j][0] = tmp_velocity[0];
				velocity[j][1] = tmp_velocity[1];
				value[j] = tmp_value;
				pBestPos[j][0] = tmp_pBestPos[0];
				pBestPos[j][1] = tmp_pBestPos[1];
				pBestValue[j] = tmp_pBestValue;	
			}
		}
	}
}

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

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

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

// 最適化結果のテキスト表示
void print_particles()
{
	for (int i = 0; i < SWARM_SIZE; i++) {
		double xx = get_particle_x(i);
		double yy = get_particle_y(i);
		double zz = get_particle_z(i);
		cout << "x:" << xx << " y:" << yy << " z:" << zz << endl;
	}
}

void create_optimize_data()
{
	for (int i = 0; i < ITERATION; i++) {
		cout << endl << i << endl;
		if (i > 0)
			move_particles();
		print_particles();
		for (int j = 0; j < SWARM_SIZE; j++) {
			x[i][j] = get_particle_x(j);
			y[i][j] = get_particle_y(j);
			z[i][j] = get_particle_z(j);
		}
	}

	sort_particles(); // 結果のソート
	cout << endl << "sorted result" << endl;
	print_particles(); // ソート結果の出力
}

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

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

	fout << "ITERATION " << ITERATION << endl;
	cout << "ITERATION " << ITERATION << 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++) {
		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();
}

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

	set_function_conditions();

	create_particles(); // 粒子群の生成

	create_optimize_data();

	delete_particles();

	write_data();
}

five_well_potential_ps.cpp

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

// 0以上1以下の実数乱数
#define RAND_01 ((double)rand() / RAND_MAX)

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

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 INERTIA = 0.75; // 慣性係数
double ACCEL_P = 1.0; // 加速係数（パーソナルベスト）
double ACCEL_G = 1.0; // 加速係数（グローバルベスト）

double** pos = nullptr; // 位置
double** velocity = nullptr; // 速度
double* value = nullptr; // 評価値
double** pBestPos = nullptr; // パーソナルベストの位置
double* pBestValue = nullptr; // パーソナルベストの評価値

double gBestPos[2]; // グローバルベストの位置
double gBestValue; // グローバルベストの評価値

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

// 評価値を算出し，パーソナルベストを更新する
void evaluate(int target)
{
	// 評価値を算出する
	value[target] = function(pos[target][0], pos[target][1]);

	// パーソナルベストを更新する
	if (pBestValue[target] > value[target]) {
		for (int i = 0; i < 2; i++)
			pBestPos[target][i] = pos[target][i];
		pBestValue[target] = value[target];
	}
}

void create_particle()
{
	pos = new double*[SWARM_SIZE];
	for (int i = 0; i < SWARM_SIZE; i++) {
		pos[i] = new double[2];
		pos[i][0] = X_MIN + (X_MAX - X_MIN) * RAND_01;
		pos[i][1] = Y_MIN + (Y_MAX - Y_MIN) * RAND_01;
	}

	velocity = new double*[SWARM_SIZE];
	for (int i = 0; i < SWARM_SIZE; i++) {
		velocity[i] = new double[2];
		velocity[i][0] = X_MIN + (X_MAX - X_MIN) * RAND_01;
		velocity[i][1] = Y_MIN + (Y_MAX - Y_MIN) * RAND_01;
	}

	value = new double[SWARM_SIZE];

	pBestPos = new double*[SWARM_SIZE];
	for (int i = 0; i < SWARM_SIZE; i++)
		pBestPos[i] = new double[2];

	pBestValue = new double[SWARM_SIZE];
	for (int i = 0; i < SWARM_SIZE; i++)
		pBestValue[i] = DBL_MAX;

	for (int i = 0; i < SWARM_SIZE; i++)
		evaluate(i);
}

// 個々の粒子を移動する
void move_particle(int target)
{
	for (int j = 0; j < 2; j++) {
		velocity[target][j] = INERTIA * velocity[target][j] 
			+ ACCEL_G * (gBestPos[j] - pos[target][j]) * RAND_01
			+ ACCEL_P * (pBestPos[target][j] - pos[target][j]) * RAND_01;
		pos[target][j] += velocity[target][j];
	}
	pos[target][0] = max1(min1(pos[target][0], X_MAX), X_MIN);
	pos[target][1] = max1(min1(pos[target][1], Y_MAX), Y_MIN);
	evaluate(target);
}

// 粒子群の生成
void create_particles()
{
	create_particle();

	int best = 0;
	for (int i = 0; i < SWARM_SIZE; i++)
		if (value[best] > value[i])
			best = i;

	for (int i = 0; i < 2; i++)
		gBestPos[i] = pos[best][i];

	gBestValue = value[best];
}

// 粒子群を初期化するため、粒子群を削除する
void delete_particles()
{
	for (int i = 0; i < SWARM_SIZE; i++)
		delete pos[i];
	delete []pos;
	pos = nullptr;

	for (int i = 0; i < SWARM_SIZE; i++)
		delete velocity[i];
	delete []velocity;
	velocity = nullptr;

	delete []value;
	value = nullptr;

	for (int i = 0; i < SWARM_SIZE; i++)
		delete pBestPos[i];
	delete []pBestPos;
	pBestPos = nullptr;

	delete []pBestValue;
	pBestValue = nullptr;
}

// 粒子を移動し，グローバルベストを更新する
void move_particles()
{
	// すべての粒子を移動する
	int best = -1;
	for (int j = 0; j < SWARM_SIZE; j++) {
		move_particle(j);
		if (gBestValue > value[j])
			best = j;
	}

	// グローバルベストを更新する
	if (best != -1) {
		gBestPos[0] = pos[best][0];
		gBestPos[1] = pos[best][1];
		gBestValue = value[best];
	}
}

// パーティクルを評価値の高い順にソートする
void sort_particles()
{
	for (int i = 0; i < SWARM_SIZE - 1; i++) {
		for (int j = i + 1; j < SWARM_SIZE; j++) {
			if (value[i] > value[j]) {
				double tmp_pos[2];
				double tmp_velocity[2];
				double tmp_value;
				double tmp_pBestPos[2];
				double tmp_pBestValue;

				tmp_pos[0] = pos[i][0];
				tmp_pos[1] = pos[i][1];
				tmp_velocity[0] = velocity[i][0];
				tmp_velocity[1] = velocity[i][1];
				tmp_value = value[i];
				tmp_pBestPos[0] = pBestPos[i][0];
				tmp_pBestPos[1] = pBestPos[i][1];
				tmp_pBestValue = pBestValue[i];

				pos[i][0] = pos[j][0];
				pos[i][1] = pos[j][1];
				velocity[i][0] = velocity[j][0];
				velocity[i][1] = velocity[j][1];
				value[i] = value[j];
				pBestPos[i][0] = pBestPos[j][0];
				pBestPos[i][1] = pBestPos[j][1];
				pBestValue[i] = pBestValue[j];

				pos[j][0] = tmp_pos[0];
				pos[j][1] = tmp_pos[1];
				velocity[j][0] = tmp_velocity[0];
				velocity[j][1] = tmp_velocity[1];
				value[j] = tmp_value;
				pBestPos[j][0] = tmp_pBestPos[0];
				pBestPos[j][1] = tmp_pBestPos[1];
				pBestValue[j] = tmp_pBestValue;	
			}
		}
	}
}

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

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

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

// 最適化結果のテキスト表示
void print_particles()
{
	for (int i = 0; i < SWARM_SIZE; i++) {
		double xx = get_particle_x(i);
		double yy = get_particle_y(i);
		double zz = get_particle_z(i);
		cout << "x:" << xx << " y:" << yy << " z:" << zz << endl;
	}
}

void create_optimize_data()
{
	for (int i = 0; i < ITERATION; i++) {
		cout << endl << i << endl;
		if (i > 0)
			move_particles();
		print_particles();
		for (int j = 0; j < SWARM_SIZE; j++) {
			x[i][j] = get_particle_x(j);
			y[i][j] = get_particle_y(j);
			z[i][j] = get_particle_z(j);
		}
	}

	sort_particles(); // 結果のソート
	cout << endl << "sorted result" << endl;
	print_particles(); // ソート結果の出力
}

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

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

	fout << "ITERATION " << ITERATION << endl;
	cout << "ITERATION " << ITERATION << 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++) {
		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();
}

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

	set_function_conditions();

	create_particles(); // 粒子群の生成

	create_optimize_data();

	delete_particles();

	write_data();
}

himmelblau_ps.cpp

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

// 0以上1以下の実数乱数
#define RAND_01 ((double)rand() / RAND_MAX)

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

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 INERTIA = 0.75; // 慣性係数
double ACCEL_P = 1.0; // 加速係数（パーソナルベスト）
double ACCEL_G = 1.0; // 加速係数（グローバルベスト）

double** pos = nullptr; // 位置
double** velocity = nullptr; // 速度
double* value = nullptr; // 評価値
double** pBestPos = nullptr; // パーソナルベストの位置
double* pBestValue = nullptr; // パーソナルベストの評価値

double gBestPos[2]; // グローバルベストの位置
double gBestValue; // グローバルベストの評価値

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

// 評価値を算出し，パーソナルベストを更新する
void evaluate(int target)
{
	// 評価値を算出する
	value[target] = function(pos[target][0], pos[target][1]);

	// パーソナルベストを更新する
	if (pBestValue[target] > value[target]) {
		for (int i = 0; i < 2; i++)
			pBestPos[target][i] = pos[target][i];
		pBestValue[target] = value[target];
	}
}

void create_particle()
{
	pos = new double*[SWARM_SIZE];
	for (int i = 0; i < SWARM_SIZE; i++) {
		pos[i] = new double[2];
		pos[i][0] = X_MIN + (X_MAX - X_MIN) * RAND_01;
		pos[i][1] = Y_MIN + (Y_MAX - Y_MIN) * RAND_01;
	}

	velocity = new double*[SWARM_SIZE];
	for (int i = 0; i < SWARM_SIZE; i++) {
		velocity[i] = new double[2];
		velocity[i][0] = X_MIN + (X_MAX - X_MIN) * RAND_01;
		velocity[i][1] = Y_MIN + (Y_MAX - Y_MIN) * RAND_01;
	}

	value = new double[SWARM_SIZE];

	pBestPos = new double*[SWARM_SIZE];
	for (int i = 0; i < SWARM_SIZE; i++)
		pBestPos[i] = new double[2];

	pBestValue = new double[SWARM_SIZE];
	for (int i = 0; i < SWARM_SIZE; i++)
		pBestValue[i] = DBL_MAX;

	for (int i = 0; i < SWARM_SIZE; i++)
		evaluate(i);
}

// 個々の粒子を移動する
void move_particle(int target)
{
	for (int j = 0; j < 2; j++) {
		velocity[target][j] = INERTIA * velocity[target][j] 
			+ ACCEL_G * (gBestPos[j] - pos[target][j]) * RAND_01
			+ ACCEL_P * (pBestPos[target][j] - pos[target][j]) * RAND_01;
		pos[target][j] += velocity[target][j];
	}
	pos[target][0] = max1(min1(pos[target][0], X_MAX), X_MIN);
	pos[target][1] = max1(min1(pos[target][1], Y_MAX), Y_MIN);
	evaluate(target);
}

// 粒子群の生成
void create_particles()
{
	create_particle();

	int best = 0;
	for (int i = 0; i < SWARM_SIZE; i++)
		if (value[best] > value[i])
			best = i;

	for (int i = 0; i < 2; i++)
		gBestPos[i] = pos[best][i];

	gBestValue = value[best];
}

// 粒子群を初期化するため、粒子群を削除する
void delete_particles()
{
	for (int i = 0; i < SWARM_SIZE; i++)
		delete pos[i];
	delete []pos;
	pos = nullptr;

	for (int i = 0; i < SWARM_SIZE; i++)
		delete velocity[i];
	delete []velocity;
	velocity = nullptr;

	delete []value;
	value = nullptr;

	for (int i = 0; i < SWARM_SIZE; i++)
		delete pBestPos[i];
	delete []pBestPos;
	pBestPos = nullptr;

	delete []pBestValue;
	pBestValue = nullptr;
}

// 粒子を移動し，グローバルベストを更新する
void move_particles()
{
	// すべての粒子を移動する
	int best = -1;
	for (int j = 0; j < SWARM_SIZE; j++) {
		move_particle(j);
		if (gBestValue > value[j])
			best = j;
	}

	// グローバルベストを更新する
	if (best != -1) {
		gBestPos[0] = pos[best][0];
		gBestPos[1] = pos[best][1];
		gBestValue = value[best];
	}
}

// パーティクルを評価値の高い順にソートする
void sort_particles()
{
	for (int i = 0; i < SWARM_SIZE - 1; i++) {
		for (int j = i + 1; j < SWARM_SIZE; j++) {
			if (value[i] > value[j]) {
				double tmp_pos[2];
				double tmp_velocity[2];
				double tmp_value;
				double tmp_pBestPos[2];
				double tmp_pBestValue;

				tmp_pos[0] = pos[i][0];
				tmp_pos[1] = pos[i][1];
				tmp_velocity[0] = velocity[i][0];
				tmp_velocity[1] = velocity[i][1];
				tmp_value = value[i];
				tmp_pBestPos[0] = pBestPos[i][0];
				tmp_pBestPos[1] = pBestPos[i][1];
				tmp_pBestValue = pBestValue[i];

				pos[i][0] = pos[j][0];
				pos[i][1] = pos[j][1];
				velocity[i][0] = velocity[j][0];
				velocity[i][1] = velocity[j][1];
				value[i] = value[j];
				pBestPos[i][0] = pBestPos[j][0];
				pBestPos[i][1] = pBestPos[j][1];
				pBestValue[i] = pBestValue[j];

				pos[j][0] = tmp_pos[0];
				pos[j][1] = tmp_pos[1];
				velocity[j][0] = tmp_velocity[0];
				velocity[j][1] = tmp_velocity[1];
				value[j] = tmp_value;
				pBestPos[j][0] = tmp_pBestPos[0];
				pBestPos[j][1] = tmp_pBestPos[1];
				pBestValue[j] = tmp_pBestValue;	
			}
		}
	}
}

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

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

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

// 最適化結果のテキスト表示
void print_particles()
{
	for (int i = 0; i < SWARM_SIZE; i++) {
		double xx = get_particle_x(i);
		double yy = get_particle_y(i);
		double zz = get_particle_z(i);
		cout << "x:" << xx << " y:" << yy << " z:" << zz << endl;
	}
}

void create_optimize_data()
{
	for (int i = 0; i < ITERATION; i++) {
		cout << endl << i << endl;
		if (i > 0)
			move_particles();
		print_particles();
		for (int j = 0; j < SWARM_SIZE; j++) {
			x[i][j] = get_particle_x(j);
			y[i][j] = get_particle_y(j);
			z[i][j] = get_particle_z(j);
		}
	}

	sort_particles(); // 結果のソート
	cout << endl << "sorted result" << endl;
	print_particles(); // ソート結果の出力
}

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

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

	fout << "ITERATION " << ITERATION << endl;
	cout << "ITERATION " << ITERATION << 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++) {
		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();
}

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

	set_function_conditions();

	create_particles(); // 粒子群の生成

	create_optimize_data();

	delete_particles();

	write_data();
}

rosenbrock_ps.cpp

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

// 0以上1以下の実数乱数
#define RAND_01 ((double)rand() / RAND_MAX)

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

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 INERTIA = 0.75; // 慣性係数
double ACCEL_P = 1.0; // 加速係数（パーソナルベスト）
double ACCEL_G = 1.0; // 加速係数（グローバルベスト）

double** pos = nullptr; // 位置
double** velocity = nullptr; // 速度
double* value = nullptr; // 評価値
double** pBestPos = nullptr; // パーソナルベストの位置
double* pBestValue = nullptr; // パーソナルベストの評価値

double gBestPos[2]; // グローバルベストの位置
double gBestValue; // グローバルベストの評価値

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

// 評価値を算出し，パーソナルベストを更新する
void evaluate(int target)
{
	// 評価値を算出する
	value[target] = function(pos[target][0], pos[target][1]);

	// パーソナルベストを更新する
	if (pBestValue[target] > value[target]) {
		for (int i = 0; i < 2; i++)
			pBestPos[target][i] = pos[target][i];
		pBestValue[target] = value[target];
	}
}

void create_particle()
{
	pos = new double*[SWARM_SIZE];
	for (int i = 0; i < SWARM_SIZE; i++) {
		pos[i] = new double[2];
		pos[i][0] = X_MIN + (X_MAX - X_MIN) * RAND_01;
		pos[i][1] = Y_MIN + (Y_MAX - Y_MIN) * RAND_01;
	}

	velocity = new double*[SWARM_SIZE];
	for (int i = 0; i < SWARM_SIZE; i++) {
		velocity[i] = new double[2];
		velocity[i][0] = X_MIN + (X_MAX - X_MIN) * RAND_01;
		velocity[i][1] = Y_MIN + (Y_MAX - Y_MIN) * RAND_01;
	}

	value = new double[SWARM_SIZE];

	pBestPos = new double*[SWARM_SIZE];
	for (int i = 0; i < SWARM_SIZE; i++)
		pBestPos[i] = new double[2];

	pBestValue = new double[SWARM_SIZE];
	for (int i = 0; i < SWARM_SIZE; i++)
		pBestValue[i] = DBL_MAX;

	for (int i = 0; i < SWARM_SIZE; i++)
		evaluate(i);
}

// 個々の粒子を移動する
void move_particle(int target)
{
	for (int j = 0; j < 2; j++) {
		velocity[target][j] = INERTIA * velocity[target][j] 
			+ ACCEL_G * (gBestPos[j] - pos[target][j]) * RAND_01
			+ ACCEL_P * (pBestPos[target][j] - pos[target][j]) * RAND_01;
		pos[target][j] += velocity[target][j];
	}
	pos[target][0] = max1(min1(pos[target][0], X_MAX), X_MIN);
	pos[target][1] = max1(min1(pos[target][1], Y_MAX), Y_MIN);
	evaluate(target);
}

// 粒子群の生成
void create_particles()
{
	create_particle();

	int best = 0;
	for (int i = 0; i < SWARM_SIZE; i++)
		if (value[best] > value[i])
			best = i;

	for (int i = 0; i < 2; i++)
		gBestPos[i] = pos[best][i];

	gBestValue = value[best];
}

// 粒子群を初期化するため、粒子群を削除する
void delete_particles()
{
	for (int i = 0; i < SWARM_SIZE; i++)
		delete pos[i];
	delete []pos;
	pos = nullptr;

	for (int i = 0; i < SWARM_SIZE; i++)
		delete velocity[i];
	delete []velocity;
	velocity = nullptr;

	delete []value;
	value = nullptr;

	for (int i = 0; i < SWARM_SIZE; i++)
		delete pBestPos[i];
	delete []pBestPos;
	pBestPos = nullptr;

	delete []pBestValue;
	pBestValue = nullptr;
}

// 粒子を移動し，グローバルベストを更新する
void move_particles()
{
	// すべての粒子を移動する
	int best = -1;
	for (int j = 0; j < SWARM_SIZE; j++) {
		move_particle(j);
		if (gBestValue > value[j])
			best = j;
	}

	// グローバルベストを更新する
	if (best != -1) {
		gBestPos[0] = pos[best][0];
		gBestPos[1] = pos[best][1];
		gBestValue = value[best];
	}
}

// パーティクルを評価値の高い順にソートする
void sort_particles()
{
	for (int i = 0; i < SWARM_SIZE - 1; i++) {
		for (int j = i + 1; j < SWARM_SIZE; j++) {
			if (value[i] > value[j]) {
				double tmp_pos[2];
				double tmp_velocity[2];
				double tmp_value;
				double tmp_pBestPos[2];
				double tmp_pBestValue;

				tmp_pos[0] = pos[i][0];
				tmp_pos[1] = pos[i][1];
				tmp_velocity[0] = velocity[i][0];
				tmp_velocity[1] = velocity[i][1];
				tmp_value = value[i];
				tmp_pBestPos[0] = pBestPos[i][0];
				tmp_pBestPos[1] = pBestPos[i][1];
				tmp_pBestValue = pBestValue[i];

				pos[i][0] = pos[j][0];
				pos[i][1] = pos[j][1];
				velocity[i][0] = velocity[j][0];
				velocity[i][1] = velocity[j][1];
				value[i] = value[j];
				pBestPos[i][0] = pBestPos[j][0];
				pBestPos[i][1] = pBestPos[j][1];
				pBestValue[i] = pBestValue[j];

				pos[j][0] = tmp_pos[0];
				pos[j][1] = tmp_pos[1];
				velocity[j][0] = tmp_velocity[0];
				velocity[j][1] = tmp_velocity[1];
				value[j] = tmp_value;
				pBestPos[j][0] = tmp_pBestPos[0];
				pBestPos[j][1] = tmp_pBestPos[1];
				pBestValue[j] = tmp_pBestValue;	
			}
		}
	}
}

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

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

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

// 最適化結果のテキスト表示
void print_particles()
{
	for (int i = 0; i < SWARM_SIZE; i++) {
		double xx = get_particle_x(i);
		double yy = get_particle_y(i);
		double zz = get_particle_z(i);
		cout << "x:" << xx << " y:" << yy << " z:" << zz << endl;
	}
}

void create_optimize_data()
{
	for (int i = 0; i < ITERATION; i++) {
		cout << endl << i << endl;
		if (i > 0)
			move_particles();
		print_particles();
		for (int j = 0; j < SWARM_SIZE; j++) {
			x[i][j] = get_particle_x(j);
			y[i][j] = get_particle_y(j);
			z[i][j] = get_particle_z(j);
		}
	}

	sort_particles(); // 結果のソート
	cout << endl << "sorted result" << endl;
	print_particles(); // ソート結果の出力
}

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

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

	fout << "ITERATION " << ITERATION << endl;
	cout << "ITERATION " << ITERATION << 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++) {
		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();
}

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

	set_function_conditions();

	create_particles(); // 粒子群の生成

	create_optimize_data();

	delete_particles();

	write_data();
}

rastrigin_ps.cpp

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

// 0以上1以下の実数乱数
#define RAND_01 ((double)rand() / RAND_MAX)

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

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 INERTIA = 0.75; // 慣性係数
double ACCEL_P = 1.0; // 加速係数（パーソナルベスト）
double ACCEL_G = 1.0; // 加速係数（グローバルベスト）

double** pos = nullptr; // 位置
double** velocity = nullptr; // 速度
double* value = nullptr; // 評価値
double** pBestPos = nullptr; // パーソナルベストの位置
double* pBestValue = nullptr; // パーソナルベストの評価値

double gBestPos[2]; // グローバルベストの位置
double gBestValue; // グローバルベストの評価値

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

// 評価値を算出し，パーソナルベストを更新する
void evaluate(int target)
{
	// 評価値を算出する
	value[target] = function(pos[target][0], pos[target][1]);

	// パーソナルベストを更新する
	if (pBestValue[target] > value[target]) {
		for (int i = 0; i < 2; i++)
			pBestPos[target][i] = pos[target][i];
		pBestValue[target] = value[target];
	}
}

void create_particle()
{
	pos = new double*[SWARM_SIZE];
	for (int i = 0; i < SWARM_SIZE; i++) {
		pos[i] = new double[2];
		pos[i][0] = X_MIN + (X_MAX - X_MIN) * RAND_01;
		pos[i][1] = Y_MIN + (Y_MAX - Y_MIN) * RAND_01;
	}

	velocity = new double*[SWARM_SIZE];
	for (int i = 0; i < SWARM_SIZE; i++) {
		velocity[i] = new double[2];
		velocity[i][0] = X_MIN + (X_MAX - X_MIN) * RAND_01;
		velocity[i][1] = Y_MIN + (Y_MAX - Y_MIN) * RAND_01;
	}

	value = new double[SWARM_SIZE];

	pBestPos = new double*[SWARM_SIZE];
	for (int i = 0; i < SWARM_SIZE; i++)
		pBestPos[i] = new double[2];

	pBestValue = new double[SWARM_SIZE];
	for (int i = 0; i < SWARM_SIZE; i++)
		pBestValue[i] = DBL_MAX;

	for (int i = 0; i < SWARM_SIZE; i++)
		evaluate(i);
}

// 個々の粒子を移動する
void move_particle(int target)
{
	for (int j = 0; j < 2; j++) {
		velocity[target][j] = INERTIA * velocity[target][j] 
			+ ACCEL_G * (gBestPos[j] - pos[target][j]) * RAND_01
			+ ACCEL_P * (pBestPos[target][j] - pos[target][j]) * RAND_01;
		pos[target][j] += velocity[target][j];
	}
	pos[target][0] = max1(min1(pos[target][0], X_MAX), X_MIN);
	pos[target][1] = max1(min1(pos[target][1], Y_MAX), Y_MIN);
	evaluate(target);
}

// 粒子群の生成
void create_particles()
{
	create_particle();

	int best = 0;
	for (int i = 0; i < SWARM_SIZE; i++)
		if (value[best] > value[i])
			best = i;

	for (int i = 0; i < 2; i++)
		gBestPos[i] = pos[best][i];

	gBestValue = value[best];
}

// 粒子群を初期化するため、粒子群を削除する
void delete_particles()
{
	for (int i = 0; i < SWARM_SIZE; i++)
		delete pos[i];
	delete []pos;
	pos = nullptr;

	for (int i = 0; i < SWARM_SIZE; i++)
		delete velocity[i];
	delete []velocity;
	velocity = nullptr;

	delete []value;
	value = nullptr;

	for (int i = 0; i < SWARM_SIZE; i++)
		delete pBestPos[i];
	delete []pBestPos;
	pBestPos = nullptr;

	delete []pBestValue;
	pBestValue = nullptr;
}

// 粒子を移動し，グローバルベストを更新する
void move_particles()
{
	// すべての粒子を移動する
	int best = -1;
	for (int j = 0; j < SWARM_SIZE; j++) {
		move_particle(j);
		if (gBestValue > value[j])
			best = j;
	}

	// グローバルベストを更新する
	if (best != -1) {
		gBestPos[0] = pos[best][0];
		gBestPos[1] = pos[best][1];
		gBestValue = value[best];
	}
}

// パーティクルを評価値の高い順にソートする
void sort_particles()
{
	for (int i = 0; i < SWARM_SIZE - 1; i++) {
		for (int j = i + 1; j < SWARM_SIZE; j++) {
			if (value[i] > value[j]) {
				double tmp_pos[2];
				double tmp_velocity[2];
				double tmp_value;
				double tmp_pBestPos[2];
				double tmp_pBestValue;

				tmp_pos[0] = pos[i][0];
				tmp_pos[1] = pos[i][1];
				tmp_velocity[0] = velocity[i][0];
				tmp_velocity[1] = velocity[i][1];
				tmp_value = value[i];
				tmp_pBestPos[0] = pBestPos[i][0];
				tmp_pBestPos[1] = pBestPos[i][1];
				tmp_pBestValue = pBestValue[i];

				pos[i][0] = pos[j][0];
				pos[i][1] = pos[j][1];
				velocity[i][0] = velocity[j][0];
				velocity[i][1] = velocity[j][1];
				value[i] = value[j];
				pBestPos[i][0] = pBestPos[j][0];
				pBestPos[i][1] = pBestPos[j][1];
				pBestValue[i] = pBestValue[j];

				pos[j][0] = tmp_pos[0];
				pos[j][1] = tmp_pos[1];
				velocity[j][0] = tmp_velocity[0];
				velocity[j][1] = tmp_velocity[1];
				value[j] = tmp_value;
				pBestPos[j][0] = tmp_pBestPos[0];
				pBestPos[j][1] = tmp_pBestPos[1];
				pBestValue[j] = tmp_pBestValue;	
			}
		}
	}
}

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

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

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

// 最適化結果のテキスト表示
void print_particles()
{
	for (int i = 0; i < SWARM_SIZE; i++) {
		double xx = get_particle_x(i);
		double yy = get_particle_y(i);
		double zz = get_particle_z(i);
		cout << "x:" << xx << " y:" << yy << " z:" << zz << endl;
	}
}

void create_optimize_data()
{
	for (int i = 0; i < ITERATION; i++) {
		cout << endl << i << endl;
		if (i > 0)
			move_particles();
		print_particles();
		for (int j = 0; j < SWARM_SIZE; j++) {
			x[i][j] = get_particle_x(j);
			y[i][j] = get_particle_y(j);
			z[i][j] = get_particle_z(j);
		}
	}

	sort_particles(); // 結果のソート
	cout << endl << "sorted result" << endl;
	print_particles(); // ソート結果の出力
}

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

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

	fout << "ITERATION " << ITERATION << endl;
	cout << "ITERATION " << ITERATION << 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++) {
		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();
}

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

	set_function_conditions();

	create_particles(); // 粒子群の生成

	create_optimize_data();

	delete_particles();

	write_data();
}

出力結果 (sphere_ps)

x * x + y * y
ITERATION 100
SWARM_SIZE 20
X_MAX 10
X_MIN -10
Y_MAX 10
Y_MIN -10
Z_MAX 200
Z_MIN 0
DATA
-3.41533 -1.8656 15.1449
-8.44966 -5.36607 100.191
-7.20389 -6.16382 89.8888
0.848109 9.2938 87.094
5.89465 9.13694 118.231
-6.29749 6.84194 86.4706
7.76299 -8.56136 133.561
8.73714 -1.07517 77.4937
-9.06735 -9.26267 168.014
3.20231 2.07801 14.5729
-9.07163 -4.91928 106.494
-0.902432 1.87597 4.33366
6.30787 -7.05618 89.579
8.54305 -0.685751 73.4539
-7.97601 -4.35163 82.5535
-0.965911 -2.46193 6.99407
-8.26472 2.93435 76.916
9.20835 -2.94656 93.4759
4.67879 4.83322 45.2511
-4.26008 -2.46803 24.2395

-6.89789 0.407926 47.7473
-5.05205 -7.03487 75.0126
-1.4152 3.37662 13.4044
6.09241 10 137.117
-0.803967 1.28332 2.29328
-3.41946 0.568544 12.0159
2.19543 -8.9841 85.5339
-1.5685 -5.3179 30.7402
-3.5741 2.57375 19.3984
5.4449 -5.31423 57.888
-10 2.29563 105.27
-7.33993 8.03835 118.49
-1.39304 -2.3947 7.67513
10 -3.18303 110.132
-5.25095 -9.08097 110.036
-7.15512 -3.82645 65.8375
-7.34208 7.60831 111.792
8.99184 2.6261 87.7496
-2.61526 7.84344 68.3591
0.443572 -0.722695 0.719044

-2.22932 -0.750848 5.53366
-1.37577 -6.49705 44.1044
3.19799 9.31944 97.079
2.9103 4.95636 33.0353
-4.77776 -6.36821 63.3811
2.12427 -4.85213 28.0557
-3.20227 -4.55867 31.036
-8.08707 -4.48828 85.5453
2.32189 10 105.391
2.65373 -5.26762 34.7901
-10 6.54635 142.855
-6.58947 0.380322 43.5658
-5.53361 2.08563 34.9707
3.26676 -3.64069 23.9264
-1.78289 -3.45404 15.109
-8.16884 -2.7245 74.1528
-5.98323 -0.71654 36.3124
3.98334 4.07992 32.5128
-5.83397 5.95161 69.4568
3.97131 0.586307 16.1151