以下の論文のPSOアルゴリズムに従った。
In [1]:
%matplotlib inline
import numpy as np
import pylab as pl
import math
from sympy import *
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from mpl_toolkits.mplot3d import Axes3D
In [182]:
def TSP_map(N): #100×100の正方格子内にN個の点を配置する関数
TSP_map = []
X = [i for i in range(100)]
Y = [i for i in range(100)]
x = np.array([])
y = np.array([])
for i in range(N):
x = np.append(x, np.random.choice(X))
y = np.append(y, np.random.choice(Y))
for i in range(N):
TSP_map.append([x[i], y[i]])
return TSP_map
In [400]:
class PSO:
def __init__(self, N, pN, omega, alpha, beta):
self.N = N
self.pN = pN
self.omega = omega
self.alpha = alpha
self.beta = beta
self.city = TSP_map(N)
def initialize(self):
ptcl = np.array([])
for i in range(self.pN):
a = np.random.choice([i for i in range(self.N - 1)])
b = np.random.choice([i for i in range(a, self.N)])
V = [a, b]
ptcl = np.append(ptcl, particle(i, V, self.N, self.omega, self.alpha, self.beta))
self.ptcl = ptcl
return self.ptcl
def one_simulate(self):
for i in range(self.pN):
self.ptcl[i].SS_id()
self.ptcl[i].SS_gd(self.p_gd_X)
self.ptcl[i].new_V()
self.ptcl[i].new_X()
self.ptcl[i].P_id(self.city)
def simulate(self, sim_num):
for i in range(self.pN):
self.ptcl[i].initial(self.city)
self.p_gd_X = self.P_gd()
for i in range(sim_num):
self.one_simulate()
self.p_gd_X = self.P_gd()
self.p_gd_X = self.P_gd()
return self.p_gd_X
def P_gd(self):
P_gd = self.ptcl[0].p_id
no = 0
for i in range(self.pN):
if P_gd < self.ptcl[i].p_id:
P_gd = self.ptcl[i].p_id
self.no = i
return self.ptcl[self.no].p_id_X
In [417]:
class particle:
def __init__(self, No, V, num_city, omega, alpha, beta): #noは粒子の番号（ナンバー）である。
self.No = No
self.V = V
self.num_city = num_city
self.omega = omega
self.alpha = alpha
self.beta = beta
self.X = self.init_X()
def initial(self, city):
self.ss_id = []
self.ss_gd = []
self.P_id(city)
def init_X(self):
c = np.array([i for i in range(self.num_city)])
np.random.shuffle(c)
return c
def SO(self, V, P):
SO = []
for i in range(len(V)):
if V[i] != P[i]:
t = np.where(V == P[i])
t = int(t[0])
a = V[i]
b = V[t]
V[i] = b
V[t] = a
SO.append([i, t])
return SO
def SS_id(self):
self.ss_id = self.SO(self.X, self.p_id_X)
def SS_gd(self, p_gd_X):
self.ss_gd = self.SO(self.X, p_gd_X)
def select(self, V, v, p):
select_v = np.array([])
for i in range(len(V)):
x = np.random.choice([1, 0], p=[p, 1-p])
if x == 1:
select_v = np.append(select_v,V[i])
return select_v
def new_V(self):
V = np.array([])
self.select(V, self.V, self.omega)
self.select(V, self.SS_id, self.alpha)
self.select(V, self.ss_gd, self.beta)
self.V = V
return self.V
def new_X(self):
for i in range(len(self.V)):
j = self.V[i][0]
k = self.V[i][1]
a = self.X[j]
b = self.X[k]
self.X[j] = b
self.X[k] = a
return self.X
def P_id(self, city): #二都市間の距離を足し上げてP_idを求める関数。
P_id = 0
for i in range(self.num_city):
if i != self.num_city-1:
x1 = city[self.X[i]][0]
y1 = city[self.X[i]][1]
x2 = city[self.X[i+1]][0]
y2 = city[self.X[i+1]][1]
else:
x1 = city[self.X[i]][0]
y1 = city[self.X[i]][1]
x2 = city[self.X[0]][0]
y2 = city[self.X[0]][1]
a = np.array([x1, y1])
b = np.array([x2, y2])
u = b - a
p = np.linalg.norm(u)
P_id += p
if P_id < self.P_id:
self.p_id = P_id
self.p_id_X = self.X
return self.p_id
10都市で行う。パラメータは左から順に、（都市の数、粒子の数、前期の速度の影響率、localベストの影響率、globalベストの影響率)である。
In [418]:
pso = PSO(10, 10, 0.3, 0.1, 0.3)
都市の座標とプロット。
In [419]:
pso.city
Out[419]:
In [420]:
x = []
y = []
for i in range(len(pso.city)):
x.append(pso.city[i][0])
y.append(pso.city[i][1])
plt.scatter(x, y)
Out[420]:
粒子の初期化。
In [421]:
pso.initialize()
Out[421]:
100回シミュレーションした。
In [422]:
pso.simulate(100)
In [415]:
x = []
y = []
for i in pso.p_gd_X:
x.append(pso.city[i][0])
y.append(pso.city[i][1])
plt.plot(x, y)
Out[415]:
In [416]:
pso.no
Out[416]:
In [ ]:
In [ ]: