Open ntaku256 opened 9 months ago
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as anime
import openpyxl
import random
from sklearn.cluster import KMeans
from mpl_toolkits.mplot3d import Axes3D
from scipy.spatial.distance import cdist
from benchmark_function import BenchmarkFunction as BF
bf = BF()
def SelectOptimization():
optimization_problems = {
"sphere": [-5.12, 5.12],
"ellipsoid": [-5.12, 5.12],
"k_tablet": [-5.12, 5.12],
"rosenbrock_star": [-2.048, 2.048],
"rosenbrock_chain": [-2.048, 2.048],
"bohachevsky": [-5.12, 5.12],
"ackley": [-32.768, 32.768],
"schaffer": [-100, 100],
"rastrigin": [-5.12, 5.12],
}
print("次の最適化問題から選んでください:")
for index, problem in enumerate(optimization_problems.keys(), start=1):
print(f"{index}. {problem}")
while True:
try:
user_input = int(input("選択したい最適化問題の番号を入力してください: "))
if 1 <= user_input <= len(optimization_problems):
problem_name = list(optimization_problems.keys())[user_input - 1]
bounds = optimization_problems[problem_name]
print(f"選択された最適化問題: {problem_name}")
print(f"範囲: {bounds}")
return problem_name, optimization_problems[problem_name][1]
else:
print("無効な選択です。1から{}の番号を入力してください。".format(len(optimization_problems)))
except ValueError:
print("無効な入力です。番号を整数として入力してください。")
# 評価値の計算
def Evaluate(problem_name, vector):
match problem_name:
case "sphere":
return bf.sphere(vector)
case "ellipsoid":
return bf.ellipsoid(vector)
case "k_tablet":
return bf.k_tablet(vector)
case "rosenbrock_star":
return bf.rosenbrock_star(vector)
case "rosenbrock_chain":
return bf.rosenbrock_chain(vector)
case "bohachevsky":
return bf.bohachevsky(vector)
case "ackley":
return bf.ackley(vector)
case "schaffer":
return bf.schaffer(vector)
case "rastrigin":
return bf.rastrigin(vector)
# 配列の中から、値の大きさに応じた確率(大きさの比)でインデックスを選択する
def roulett(table):
rand = np.random.uniform(0.0, np.sum(table))
sum = 0
for i in range(len(table)):
sum += table[i]
if sum > rand:
return i
# 配列の中から、値の大きさに応じた確率(逆数の比)でインデックスを選択する
# ただし、0付近で値が極端に大きくなるのでオフセットを使用する
# オフセットの値は、求める解の精度と扱える精度を考慮して決める
def roulettMin(table, offset=1e-300):
total = 0
for i in range(len(table)):
total += 1 / (offset + table[i])
rand = random.uniform(0.0, total)
sum = 0
for i in range(len(table)):
sum += 1 / (offset + table[i])
if sum > rand:
return i
# 粒子の位置が範囲外に出た場合は、範囲内の最大値にする
def filtIndivisual(vector, r):
for i in range(len(vector)):
vector[i] = max(-r, min(r, vector[i]))
# 粒子の位置を初期化
def InitIndivisual(r):
return np.random.uniform(-r, r, (1 * 2))
class Instabae:
n_iters = None # 試行回数
n_indivisuals = None # 粒子の数
n_clusters = None # クラスタ数
r = None # 探索範囲
w = None
particles_position = None # 粒子の位置
particles_velocity = None # 粒子の速度
particles_value = None # 粒子の評価値
particles_label = None
particles_color = None
particles_best_positon = None # 粒子内の最適解の位置
particles_best_value = None # 粒子内の最適解の評価値
cluster_center_position = None
cluster_center_value = None
problem_name = None
top_percent = 0.05
max_clusters = 20
def __init__(self, problem_name, n_iters, n_indivisuals, r, w):
self.problem_name = problem_name
self.n_iters = n_iters
self.n_indivisuals = n_indivisuals
self.r = r
self.w = w
self.InitIndivisuals()
self.Evaluate()
def InitIndivisuals(self):
self.particles_position = np.array([InitIndivisual(self.r) for _ in range(self.n_indivisuals)])
self.particles_velocity = np.zeros_like(self.particles_position)
self.particles_value = np.array([float("inf") for _ in range(self.n_indivisuals)])
self.particles_best_positon = np.zeros_like(self.particles_position)
self.particles_best_value = np.array([float("inf") for _ in range(self.n_indivisuals)])
self.particles_color = np.array([(0.0, 0.0, 0.0, 0.0) for _ in range(self.n_indivisuals)])
self.cluster_center_position = np.zeros_like(self.particles_position)
self.cluster_center_value = np.array([float("inf") for _ in range(self.n_indivisuals)])
def InitClustering(self):
# 上位any%の粒子のインデックスを取得
n_top_individuals = int(self.n_indivisuals * self.top_percent)
top_indices = np.argsort(self.particles_value)[:n_top_individuals]
self.n_clusters = n_top_individuals
# 上位any%の粒子の位置を初期重心として抽出
self.cluster_center_position = self.particles_position[top_indices]
# 各粒子の重心への距離を計算し、最も近い重心に割り当て
distances = cdist(self.particles_position, self.cluster_center_position, metric="euclidean")
self.particles_label = np.argmin(distances, axis=1)
# カラーマップからクラスタ数に応じた色を生成
cmap = plt.cm.get_cmap("tab20", self.n_clusters) # 'tab20'で最大20色の見分けやすいカラーマップを取得
colors = [cmap(i) for i in range(self.n_clusters)]
for cluster_id in range(self.n_clusters):
# クラスタIDに属するすべての粒子に対して色を一括で割り当てる
self.particles_color[self.particles_label == cluster_id] = colors[cluster_id]
def Evaluate(self):
self.particles_value = Evaluate(self.problem_name, self.particles_position)
for i in range(self.n_indivisuals):
if self.particles_value[i] < self.particles_best_value[i]:
self.particles_best_value[i] = self.particles_value[i]
self.particles_best_positon[i] = self.particles_position[i]
def UpdateVectors(self):
for i in range(self.n_indivisuals):
label = self.particles_label[i]
r1 = np.random.uniform(0, 1, self.particles_position[0].shape) * 0.7
r2 = np.random.uniform(0, 1, self.particles_position[0].shape) * 0.5
center_vector = self.cluster_center_position[label] - self.particles_position[i]
local_vector = self.particles_best_positon[i] - self.particles_position[i]
self.particles_velocity[i] = self.w * self.particles_velocity[i] + r1 * center_vector + r2 * local_vector
self.particles_position[i] = self.particles_position[i] + self.particles_velocity[i]
filtIndivisual(self.particles_position[i], self.r)
for j in range(len(self.particles_position[i])):
self.particles_position[i][j] = max(-r, min(r, self.particles_position[i][j]))
def Clustering(self):
for i in range(self.n_clusters):
for j in range(self.n_indivisuals):
if self.particles_label[j] == i and self.particles_value[j] < self.cluster_center_value[i]:
self.cluster_center_value[i] = self.particles_value[j]
self.cluster_center_position[i] = self.particles_position[j]
def Run(self):
fig = plt.figure()
ax = Axes3D(fig)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("f")
fig.add_axes(ax)
X = np.linspace(-self.r, self.r, 100)
Y = np.linspace(-self.r, self.r, 100)
XX, YY = np.meshgrid(X, Y)
d = XX.shape
input_array = np.vstack((XX.flatten(), YY.flatten())).T
Z = Evaluate(self.problem_name, input_array)
ZZ = Z.reshape(d)
imgs = []
initw = self.w
self.InitClustering()
for i in range(self.n_iters):
if i < 10 or (i + 1) % 50 == 0:
title = ax.text(0, 0, 1.61 * max(Z), "%s" % (str(i + 1)), size=20, zorder=1, color="k")
scatter_vector = ax.scatter(self.particles_position.T[0], self.particles_position.T[1], self.particles_value, c=self.particles_color, s=4, alpha=1)
scatter_center_vector = ax.scatter(self.cluster_center_position.T[0], self.cluster_center_position.T[1], Evaluate(self.problem_name, self.cluster_center_position), c="black", s=8, alpha=1)
scatter_func = ax.plot_surface(XX, YY, ZZ, rstride=1, cstride=1, cmap=plt.cm.coolwarm, alpha=0.3)
imgs.append([scatter_func, scatter_center_vector, scatter_vector, title])
self.UpdateVectors()
self.Evaluate()
self.Clustering()
# self.w = self.w - ((initw - 0.4) / self.n_iters)
print(min(self.particles_best_value))
ani = anime.ArtistAnimation(fig, imgs, interval=400)
ani.save("benchmark_function/GIF/pso/pso.gif", writer="imagemagick")
plt.show()
# 2評価
# 3クラスタリング
# 4粒子の行動
return
if __name__ == "__main__":
n_iters = 1000
n_indivisuals = 150
w = 0.9
problem_name, r = SelectOptimization()
instabae = Instabae(problem_name, n_iters, n_indivisuals, r, w)
instabae.Run()
https://1.gigafile.nu/0828-g617722f93590095479c863745dc0aa6f