2025.7.29 Simulated Annealing法による1D-関数の最適化

最小化を行う1変数関数：

$ y(x) = \sin(x+1)+\sin(2x+2)+\sin(4x+3)+\sin(8x)

簡潔版から示す。

code:p.py

import numpy as np

import matplotlib.pyplot as plt

def func(x):

return np.sin(x+1) + np.sin(2*x+2) + np.sin(4*x+3) + np.sin(8*x)

def create_x():

def solve_sa(func, rho, iter_max):

T = 1.

x_now = create_x()

f_now = func(x_now)

x_best, f_best = x_now, f_now

iter_count = 0

while True:

T *= rho

x_new = create_x()

f_new = func(x_new)

metro_rand = np.random.rand()

metro_exp = np.exp(- (np.fabs(f_new - f_now) / T))

if (f_new < f_now) or (metro_rand < metro_exp):

x_now, f_now = x_new, f_new

if f_new < f_best:

x_best, f_best = x_new, f_new

if iter_count == iter_max:

print('計算回数が上限に達したので計算を打ち切った', T, iter_count)

break

if T <= eps:

print('Tが小さくなったので計算を打ち切った', T, iter_count)

break

iter_count += 1

return x_best, f_best

############################

eps = 10e-10 # Tの最小値

rho = 0.9 # Tの減少率

iter_max = 10000 # 繰り返し計算の上限

x_range = 0, 12 # 探索幅

x_best, f_best = solve_sa(func, rho, iter_max)

print('x best', x_best)

print('f best', f_best)

'''

Tが小さくなったので計算を打ち切った 9.677749120240556e-10 196

x best 3.7277164077321214

f best -2.837043805259503

'''

上例のように入力データの値が都合良く0~2,3 程度のオーダーで有ることは少ない。例えば画像データの場合は0-255で、扱われる数値の総量は画素数によって定まる。

温度Tの範囲は(0, 1)として、\Delta E の値を正規化することにより、パラメータ調整の手間を減らしてみた。

code:sa2.py

import numpy as np

import matplotlib.pyplot as plt

def func(x):

return ((np.sin(x+1) + np.sin(2*x+2) + np.sin(4*x+3) + np.sin(8*x))+2)*0.001

def create_x():

def solve_sa(func, rho, iter_max):

global x_hist, y_hist, T_hist, chk_hist, x_escape_hist, y_escape_hist

T = 1.

# x_now = create_x()

x_now = 6 # xの初期値

y_now = func(x_now)

x_best, y_best = x_now, y_now

iter_count = 0

x_hist.append(x_now)

y_hist.append(y_now)

while True:

T *= rho

x_new = create_x()

y_new = func(x_new)

metro_rand = np.random.rand()

metro_exp = np.exp(- (np.fabs((y_new - y_now)/y_now) / T))

if (y_new < y_now) or (metro_rand < metro_exp):

if y_new >= y_now:

print('escape at {} | {:.4e}({:8.4f}) -> ({:.4e}({:8.4f})'.format(iter_count, y_now, x_now, y_new, x_new))

x_escape_hist.append(x_now)

y_escape_hist.append(y_now)

x_now, y_now = x_new, y_new

if y_new < y_best:

x_best, y_best = x_new, y_new

if iter_count == iter_max:

print('計算回数が上限に達したので計算を打ち切った', T, iter_count)

break

if T <= eps:

print('Tが小さくなったので計算を打ち切った', T, iter_count)

break

x_hist.append(x_now)

y_hist.append(y_now)

iter_count += 1

T_hist.append(T)

chk_hist.append(metro_exp)

return x_best, y_best

############################

rho = 0.95

eps = 10e-10

iter_max = 1000

x_range = 0, 12

x_hist, y_hist = [], []

x_best_hist, y_best_hist = [], []

T_hist, chk_hist = [], []

x_escape_hist, y_escape_hist = [], []

x_best, y_best = solve_sa(func, rho, iter_max)

print('x best', x_best)

print('f best', y_best)

x_step = 10000

yy = func(xx)

fig = plt.figure()

fig.add_subplot(2,1,1)

plt.plot(xx, yy)

plt.plot(x_escape_hist, y_escape_hist, '*', color='red', label='espace')

plt.plot(x_hist, y_hist, '.-', label='$x$')

plt.plot(x_best, y_best, '*', label='solution')

plt.legend()

fig.add_subplot(2,1,2)

plt.plot(T_hist, label='$T$')

plt.plot(chk_hist, label='$\\exp (\\Delta)/T$')

plt.legend()

plt.show()

'''

escape at 2 | 1.1992e-03( 3.4572) -> (2.1462e-03( 7.0407)

escape at 5 | 1.4986e-03( 10.5660) -> (1.6232e-03( 8.1443)

escape at 17 | 9.9307e-05( 4.5690) -> (1.4331e-04( 10.8754)

Tが小さくなったので計算を打ち切った 9.507364444573407e-10 404

x best 9.994602764648288

f best -0.000829818296678047

'''

https://scrapbox.io/files/6889a7bfb39093532c341bc6.png