局所探索アルゴリズム・シミュレーテッドアニーリング

Image from Gyazo

ノートブックの作成

Jupyter Notebook を起動し,新規にノートブックを作成してください. ノートブックのタイトルは AI-9 とします. ノートブックの作成方法に関しては第1回の資料を参考にしてください.

最適化問題

今回も 最適化問題 を取り上げます.

最適化問題

下記の関数の値が最大となる$x$を求めよ. ただし,$x$の範囲は$0 \leq x \leq 20$とする. $$ f(x) = | \cos(x) \cdot x| $$ Image from Gyazo

最適化問題の状態空間モデルを確認しておきます.

状態空間モデル

状態は$x(0 \leq x \leq 20)$で表す. ただし,$x$は連続値ではなく,$0.1$で刻んだ離散値とする. また,現在の状態$x$から,$+0.1$,また,$-0.1$だけ増減することを行動と定義する.

  • $x' = x + 0.1$
  • $x' = x - 0.1$
初期状態は$0 \leq x \leq 20$の範囲でランダムに選んだ値とし, 目的状態は$f(x)$を最大化することである.

シミュレーテッド・アニーリング(Simulated Annealing)

シミュレーテッド・アニーリングは日本語では焼きなまし法と呼ばれる手法で, 金属加工において,金属を熱した後で徐々に冷やしながら,金属の形状を安定させる方法に由来しています. 山登り法は評価の高い方向に必ず進むのに対し,シミュレーテッド・アニーリングでは, 遷移確率$p$に応じて評価が低い場合も受け入れることにあります. これにより,局所解に陥ることなく,大局解を発見できる可能性が高くなります. 一般に,遷移確率$p$は,ステップ数(時間)が経過するとともに小さな値となり, 最終的には山登り法と同様に評価が高い方向に探索を進めます.

Image from Gyazo

Image from Gyazo

ここでは,下記の式を遷移確率$p$とします. $k$は遷移確率の最大値,$\lambda$は確率の変化を決めるパラメータであり,$i$はステップ数を表しています. 例えば,$k=50$,$\lambda=0.0005$としたとき,下図のように遷移確率は変化します. 最初は$50%$の確率で評価が低い遷移を受け入れますが, 探索が進むと確率は$0%$に近づき,高い評価のみを受け入れることになります.

$$ p = k \cdot e^{- \lambda \cdot i} $$

Image from Gyazo

実装

最初にライブラリをインポートします. 探索問題の型 を表すSearchPrblemに加え, 山登り法 を表す hill_climbingシミュレーテッド・アニーリングを表すsimulated_annealingをインポートします. また,数値計算やグラフ描画のために,NumpyMatplotibrandom, mathもインポートします.

import numpy as np
import matplotlib.pyplot as plt
import random
import math

from simpleai.search import SearchProblem, hill_climbing, simulated_annealing

最初に関数$f(x)$を定義します. ここでは,$f(x) = | \cos(x) \cdot x|$とします.

def fx(x):
    y = abs(np.cos(x) * x)
    return y

関数をグラフに描画してみましょう. グラフの描画にはMatplotlibを用います.

x = np.arange(0, 20, 0.1)
y = fx(x)
plt.plot(x, y)
plt.xlabel("x")
plt.ylabel("y")

Image from Gyazo

次に 初期状態 を定義します. 初期状態は$0$に固定します.

START = 0

Image from Gyazo

インポートしたSearchProblemクラスをオーバライドして, 最適化問題を表すHillProblemを定義します. actionsresultvalueは前回の山登り法のときと同じです.

class HillProblem(SearchProblem):
    
    def actions(self, state):

        list = []

        list.append(min(state + 0.1, 20))
        list.append(max(state - 0.1, 0))      
        
        print("{:.3f} -> [{:.3f} {:.3f}]".format(state, list[0], list[1]))

        return list 
    
    def result(self, state, action):
        return action
    
    def value(self, state):
        v = fx(state)
        return v

まずは山登り法hill_climbingを用いて探索します.

problem = HillProblem(initial_state=START)
result = hill_climbing(problem)

最終状態と評価を確認します. 最終状態は$x=0.9$,また,その評価は$f(0.9)=0559$であることがわかります.

print("{:.3f}".format(result.state))
print("{:.3f}".format(result.value))
0.900
0.559

最終状態も合わせてグラフに描画します. 局所解で探索が終了していることがわかります.

x = np.arange(0, 20, 0.1)
y = fx(x)
plt.plot(x, y)
plt.plot(result.state,  result.value,  marker="o", markersize=10)
plt.xlabel("x")
plt.ylabel("y")

Image from Gyazo

次に,シミュレーテッド・アニーリングを用いて探索します. まずは,遷移確率を表すschedule関数を定義します. $k=50$,$\lambda=0.0005$に設定しています.

def schedule(iteration, k=50, lamda=0.0005, limit=10000):
    
    p = k * math.exp(-lamda * iteration)
    
    return p

遷移確率をグラフに描画してみます. ステップ数(時間)が増加するごとに,遷移確率が$0$に近くことがわかります.

i  = np.arange(0, 10000, 1)
p = list(map(lambda x: schedule(x), i))

plt.plot(i, p)
plt.xlabel("iteration")
plt.ylabel("probability")

Image from Gyazo

それでは,シミュレーテッド・アニーリングを適用します. 最大ステップ数(繰り返し回数)は$10000$に設定しています.

problem = HillProblem(initial_state=START)
result = simulated_annealing(problem, schedule=schedule, iterations_limit=10000)

最終状態と評価を確認します. 最終状態は$x=15.8$,また,その評価は$f(15.8)=15.733$であることがわかります. 実行するごとに異なる結果になる可能性があることに注意してください.

print("{:.3f}".format(result.state))
print("{:.3f}".format(result.value))
15.800
15.733

最終状態も合わせてグラフに描画します. 大局解には到達できませんでしたが, 山登り法に比べ高い評価の状態を得られたことがわかります.

x = np.arange(0, 20, 0.1)
y = fx(x)
plt.plot(x, y)
plt.plot(result.state,  result.value,  marker="o", markersize=10)
plt.xlabel("x")
plt.ylabel("y")

Image from Gyazo

参考書籍