必要なモジュールをインポートする。
In [2]:
%matplotlib inline
import numpy as np
import math
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
rosenbrock関数である。図は2変数の場合である。
In [62]:
def rosenbrock(x):
z = 0
for i in range(1, len(x)):
z += 100*pow(x[0] - x[i], 2) + pow(1 - x[i], 2)
return z
fig = pl.figure()
ax = Axes3D(fig)
X = np.arange(-2.048, 2.048, 0.1)
X, Y = np.meshgrid(X, X)
Z = example1([X, Y])
ax.plot_surface(X, Y, Z, rstride=3, cstride=3)
pl.show()
SAとPSOを用いて実験を行う。
①先にSAで検証する。以下は提出したコードと同じ、SAのコードである。
In [487]:
class SA:
#引数にはN:変数の数、x:変数の下限と上限、init_T:温度の初期値、をとる。
def __init__(self, N=50, x=[-2.048, 2.048], init_T=100):
self.x = np.array([np.random.uniform(x[0], x[1]) for i in range(N)])
self.T = init_T
self.N = N
self.index = np.array([i for i in range(N)])
#annealingで使うための関数である。
def copy_x(self, change):
x = np.zeros(self.N)
for i in range(self.N):
if i != change:
x[i] = self.x[i]
else:
a = np.random.choice([0, 1])
if a == 0:
x[change] += 0.001*np.random.uniform(0, 2)
else:
x[change] -= 0.001*np.random.uniform(0, 2)
return x
def annealing(self):
change = np.random.choice(self.index)
#配列xのうち変更するインデックスを決定する。
x = self.copy_x(change)
#x[change]を正または負の方向にわずかに動かす。
energy = self.accept(x) - self.accept(self.x)
#変更前と変更後のエネルギー差。
#エネルギーが低くなれば値を更新、そうでなくとも一定の確率(Tに依存)で更新。
if energy < 0:
self.x = x
else:
p = pow(math.e, -energy/self.T)
b = np.random.choice([0, 1], p=[p, 1-p])
if b == 0:
self.x = x
#エネルギーの冷却
def cooling(self):
self.T = self.T*0.95
#rosenbrock関数の計算。
def accept(self, x):
cal = rosenbrock(x)
return cal
#n回繰り返す。
def simulate(self, n):
self.log = np.array([])
for i in range(n):
self.annealing()
self.cooling()
self.log = np.append(self.log, self.accept(self.x))
return self.log[-1]
SAを1000回繰り返す。下はrosenbrock関数が返す値の推移である。
In [605]:
sa = SA()
sa.simulate(1000)
Out[605]:
In [606]:
plt.plot(sa.log)
Out[606]:
何回か繰り返したが大体48を返すxの組み合わせで更新が止まってしまう（以下で再検討）。
②次にPSOで検討する。以下は提出したコードと同じ、PSOのコードである。
In [587]:
#粒子を表す。
class particle:
#引数の意味はx:粒子の位置、v:粒子の速度、alpha,beta,ganmaは速度更新で重みの調整に使われる。
def __init__(self, x, v, omega, alpha, beta, ganma):
self.x = x
self.v = v
self.omega = omega
self.alpha = alpha
self.beta = beta
self.ganma = ganma
self.lbest_x = x
self.lbest = 10000
#その粒子のローカルベストを返す関数。
def local_best(self):
cal = rosenbrock(self.x)
if self.lbest > cal:
self.lbest = cal
self.lbest_x = self.x
#vの値の更新。
def renew_v(self, gbest_x):
r_1 = np.random.uniform()
r_2 = np.random.uniform()
self.v = self.ganma*(self.omega*self.v + (self.lbest_x - self.x)*self.alpha*r_1 + (gbest_x - self.x)*self.beta*r_2)
#xの値の更新。
def renew_x(self):
self.x = self.x + self.v
#PSO全体のクラス。
class PSO:
#引数の意味はN:変数の数、pN:粒子の数、limit_x:変数（粒子の位置）の下限と上限、残りはそれぞれ対応する変数の下限と上限（値は経験的に決めた）。
def __init__(self, N=50, pN=100, limit_x=[-2.048, 2.048], limit_v=[0.0, 0.01], limit_omega=[0.8, 1.0], limit_alpha=[0.0, 1.0], limit_beta=[0.0, 1.0], limit_ganma=[0.8, 1.0]):
#Alpha is Attenuation coefficient.
#Omega is Convergence factor.
#N is the number of iteration. Exit condition is given by the number of iteration.
self.N = N
self.pN = pN
self.limit_x = limit_x
self.limit_v = limit_v
self.limit_omega = limit_omega
self.limit_alpha = limit_alpha
self.limit_ganma = limit_ganma
self.limit_beta = limit_beta
self.GBest = 10000
self.GBest_x = np.array([np.random.uniform(limit_x[0], limit_x[1]) for i in range(N)])
#以下、initialize_?関数は?の初期化を行う。
def initialize_x(self):
x = np.array([np.random.uniform(self.limit_x[0], self.limit_x[1]) for i in range(self.N)])
return x
def initialize_v(self):
v = np.array([np.random.uniform(self.limit_v[0], self.limit_v[1]) for i in range(self.N)])
return v
def initialize_omega(self):
omega = np.array([np.random.uniform(self.limit_omega[0], self.limit_omega[1]) for i in range(self.N)])
return omega
def initialize_alpha(self):
alpha = np.array([np.random.uniform(self.limit_alpha[0], self.limit_alpha[1]) for i in range(self.N)])
return alpha
def initialize_beta(self):
beta = np.array([np.random.uniform(self.limit_beta[0], self.limit_beta[1]) for i in range(self.N)])
return beta
def initialize_ganma(self):
ganma = np.array([np.random.uniform(self.limit_ganma[0], self.limit_ganma[1]) for i in range(self.N)])
return ganma
#粒子をセットする。配列に粒子オブジェクトを格納する。
def set_particles(self):
particles = np.array([])
for i in range(self.pN):
particles = np.append(particles, particle(self.initialize_x(), self.initialize_v(), self.initialize_omega(), self.initialize_alpha(), self.initialize_beta(), self.initialize_ganma()))
particles[i].local_best()
if particles[i].lbest < self.GBest:
self.GBest = particles[i].lbest
self.GBest_x = particles[i].lbest_x
particles[i].renew_v(self.GBest_x)
self.particles = particles
#n回シミュレートする。GBestはグローバルベストを示す。
def simulate(self, n):
self.set_particles()
for t in range(n):
for i in range(self.pN):
self.particles[i].renew_x()
self.particles[i].local_best()
if self.particles[i].lbest < self.GBest:
self.GBest = self.particles[i].lbest
self.GBest_x = self.particles[i].lbest_x
self.particles[i].renew_v(self.GBest_x)
return self.GBest, self.GBest_x
#1回シミュレートする。
def one_simulate(self, n):
for t in range(n):
for i in range(self.pN):
self.particles[i].renew_x()
self.particles[i].local_best()
if self.particles[i].lbest < self.GBest:
self.GBest = self.particles[i].lbest
self.GBest_x = self.particles[i].lbest_x
self.particles[i].renew_v(self.GBest_x)
return self.GBest, self.GBest_x
PSOを粒子数を100にして100回繰り返す。
In [589]:
pso = PSO()
pso.simulate(100)
Out[589]:
[test1]収束の推移を視覚化する。以下はそのための関数と、その図。
In [594]:
def test1(n):
pso = PSO()
pso.set_particles()
p = []
for i in range(n):
best, best_x = pso.one_simulate(1)
p.append(best)
plt.plot(p)
平均的な推移をみるために30回繰り返す実験を10回行っている。
In [595]:
for i in range(10):
test1(30)
平均的な推移をみるために50回繰り返す実験を10回行っている。
In [596]:
for i in range(10):
test1(50)
平均的な推移をみるために100回繰り返す実験を2回行っている。
In [597]:
for i in range(2):
test1(100)
[test2]粒子の数を200に増やしてtest1を行う。以下はそのための関数と、その図。
In [598]:
def test2(n):
pso = PSO(pN=200)
pso.set_particles()
p = []
for i in range(n):
best, best_x = pso.one_simulate(1)
p.append(best)
plt.plot(p)
In [600]:
for i in range(10):
test2(30)
[test3]粒子の数と収束の関係を確認する。以下はそのための関数と、その図。図の見出しは粒子の数。
In [601]:
def test3(n, maxpN):
pso = PSO(pN=maxpN)
pso.set_particles()
p = []
for i in range(n):
best, best_x = pso.one_simulate(1)
p.append(best)
plt.plot(p, label=maxpN)
plt.legend()
In [602]:
for i in range(1, 11):
test3(50, i*10)
In [603]:
for i in range(5, 11):
test3(50, i*10)
In [191]:
for i in range(7, 11):
test3(50, i*10)
In [604]:
for i in range(11, 20):
test3(50, i*10)
③SAとPSOを比較する。
1 SAの収束を確認する。rosenbrock関数が返す数がある値になって、xの更新が行われなくなったら試行をとめるテスト。
In [649]:
def test_sa():
sa = SA()
a = sa.simulate(10)
b = sa.simulate(10)
c = sa.simulate(10)
p = []
while c != b and c != a:
a = b
b = c
c = sa.simulate(10)
p.append(c)
plt.plot(p[5:])
print p[-1]
In [650]:
for i in range(10):
test_sa()
大体、48.9ぐらいの値を返す時点で更新をやめてしまっている。
2 PSOの収束を確認する。rosenbrock関数が返す数がある値になって、xの更新が行われなくなったら試行をとめるテストをしたかったが、収束が終わらなかったので、100回試行するごとにグローバルベストを出力させる。
In [670]:
pso = PSO()
for i in range(100):
pso.simulate(100)
print pso.GBest
100回の試行を100回繰り返したが収束は続いている。
In [671]:
for i in range(100):
pso.simulate(100)
print pso.GBest
続いて100回の試行を100回繰り返した（上と合わせて計20000回）が収束は続いている。
In [672]:
for i in range(20):
pso.simulate(500)
print pso.GBest
続いて500回の試行を20回繰り返した（上と合わせて計30000回）が収束は続いている。
In [673]:
for i in range(20):
pso.simulate(500)
print pso.GBest
続いて500回の試行を20回繰り返した（上と合わせて計40000回）が収束は続いている。
キリがないので、ここでやめる。この時点でのグローバルベストを返す粒子の位置を確認する。
In [675]:
pso.GBest_x
Out[675]:
大体、すべての値が0.52の近くにいる。これはrosenbrock関数の形的にすべての値が等しいときに、局所最適解になるからだと考えられる。
ここで10000回追加で試行すると、
In [688]:
pso.simulate(10000)
Out[688]:
となった。すべての値が0.555の近くにいる。これより、おそらくすべての値が等しい状態をできるだけ崩さす大域的最適解に向かっていると考えられる。
すべてのxが1の値をとるとき、rosenbrock関数は0となることを発見している。これが唯一の解かどうかは議論しないが、最初に述べたようにすべてのxが同じ値をとることが局所最適解であり、あるxが他のxと違う値をとると局所最適解ではなくなってしまう。そのためすべての値が1以外の値をとる局所最適解に陥ってしまい、抜け出せなくなると考えられる。SAはその局所最適解にはまっている。PSOは局所最適解に陥りながらも、少しずつすべての値が1になるようにxの値を更新しているように見える。この試行を繰り返せば大域的最適解であると予想されるすべての値が1で、rosenbrock関数が0を返す状態に収束すると考えれらるが、上記で示したように、確かに収束しているものの、収束速度が遅いため最後まで収束させることは省略する。ただし、参考のために下にxの初期値を1のまわりでのみでとった場合の試行を乗せた。以上をもってレポートを締めくくる。
[おまけ:すべてのxが1の値をとるとき、rosenbrock関数は0となることの確認]
In [687]:
c = [1.0 for i in range(50)]
rosenbrock(c)
Out[687]:
PSOでxの初期値を1のまわりでのみでとる。
In [690]:
pso = PSO(limit_x=[0.9, 1.1])
for i in range(10):
pso.simulate(1000)
print pso.GBest
In [692]:
pso.GBest_x
Out[692]:
たしかに初期値を1.0付近でとると、予想した通り、すべての値が1.0になるように収束する。
In [ ]: