ABCアルゴリズムの解説と実装【連続変数最適化】

[PR]当サイトはアフィリエイト広告を利用しています。

今回は進化計算の一つであるABCアルゴリズムについて解説とその実装を書きます。

ABCアルゴリズムは連続変数最適化で有効とされている進化計算の一つです。

ABCは「Artificial Bee Colony: ABC」からきており、人工蜂コロニーアルゴリズムとも言います。

ミツバチの集蜜作業を模したアルゴリズムで、高い性能があると言われています。

以前当ブログでは同じ進化計算アルゴリズムである遺伝的アルゴリズムの解説とPythonでのプログラミングを行いました。

また遺伝的アルゴリズムの特徴と活用例について述べました。

これらの中で連続変数のSphere関数の最適化を試してみました。Sphere関数は以下のような関数です。\(n\)は変数の数です。

\(f(\boldsymbol{x})=\sum\nolimits_{i=1}^n x_i^2\)

遺伝的アルゴリズムの活用例を述べた際、局所解が複数存在する問題の大域解の探索において遺伝的アルゴリズムが有効であるという話をしました。

ただしこうした大域解の探索において、特に連続変数の場合は遺伝的アルゴリズム以上の性能があると言われるアルゴリズムはいくつかあり、その一つがABCアルゴリズムです。

今回はABCアルゴリズムを解説し、いくつかのベンチマーク問題を解いてみて、ABCアルゴリズムがざっくりわかるようになることを目指します。

進化計算に関する知識が深まれば幸いです。

国産でちょっと安く買えるPC
マウスコンピューター

ABCアルゴリズムとは
ABCアルゴリズムのプログラミング
ABCアルゴリズムの優れた点【遺伝的アルゴリズムとの比較】
- 遺伝的アルゴリズムはランダムサーチの突然変異でしか細かい値を見つけられない
- ABCアルゴリズムは局所解に陥りにくい
遺伝的アルゴリズムおよびソルバーのエボリューショナリーとの比較
パラメータの調整
- パラメータチューニングの結果
- パラメータチューニングの結果からわかること
まとめ【連続変数の最適化ならABCアルゴリズムは優れている】

ABCアルゴリズムとは

ABCアルゴリズムの概要をお話ししていきます。

ミツバチの集蜜作業を模したアルゴリズム

ミツバチは働きバチが蜜を探しに行き、良い場所が見つかったらそれを他の働きバチに伝え、そこを重点的に探します。８の字ダンスと言われたりするこの行動を最適化アルゴリズムに取り入れた方法がABCアルゴリズムです。

ABCアルゴリズムでは「働きバチ」「追従バチ」「偵察バチ」の3つの役割をする蜂がそれぞれの役目を果たすことで最適解を探索します。

ABCアルゴリズムのアルゴリズム

ABCアルゴリズムを以下に示します。やや抽象的に書いていくので、詳細は後の節をご覧ください。

Step1：初期化をする。それぞれの働きバチが目指す場所を決める
Step2：働きバチが目指す場所の近くを探索する。今より良い場所があればそこを次の目的地にする。今より良い場所でなければ、蜜を一つ取得する
Step3：追従バチが働きバチの目標地の評価値を基に、評価値が良い場所の近くを重点的に探索する
Step4：偵察バチが新しい目標地を働きバチに示す。新しい目標地はあらかじめ決めておいた蜜の総量を取得しつくした目標地の働きバチにだけ示す。
Step5：最良値を計算して保存する
Step6：終了条件を満たしていれば終了。そうでないならStep2に戻る

Step1：初期化

以下のパラメータを決めます。

働きバチの数\(N_h\)、追従バチの数\(N_f\)、蜜の上限\(h\)を決める
\(i\)番目の働きバチが目指す場所\(\boldsymbol{x}^i\)を以下の式(1)で決める
\(i\)番目の働きバチが目指す場所の蜜の量を\(h\)に設定する

(1)\(x_{j}^{i}=x_{j}^{L}+(x_{j}^{R}-x_{j}^{L})×rand[0, 1]\)

ただし\(j = 1, 2,\ldots, m\)で\(m\)は変数の数です。

\(x_{j}^{L}\)は\(j\)番目の変数の下限、\(x_{j}^{R}\)は\(j\)番目の変数の上限です。

\(rand[0, 1]\)は0以上1未満の乱数です。

Step2：働きバチフェーズ

\(i\)番目の働きバチの新しい探索点\(\boldsymbol{y}^{i}\)を以下のように決定します。

\(y^i_j = \left\{ \begin{array}{l} x_j^i + rand[-1,1] \times ( x^i_j - x^l_j ) \quad j=k \\ x^i_j \qquad (otherwise) \end{array} \right. \)

ただし\(y^i_j\)は\(\boldsymbol{y}^{i}\)の\(j\)番目の成分です。\(k\)、\(l\)はランダムな整数で、\(k\)は1以上\(m\)以下の乱数、\(l\)は1以上\(N_h\)以下の乱数です。

各働きバチに対して一回だけ\(k\)、\(l\)を発生させます。

そして\(\boldsymbol{y}^{i}\)での目的関数の値が改善されていれば、\(\boldsymbol{y}^{i}\)を新しい目標地に設定します。新しい目標地の蜜の量を\(h\)に戻します。

改善されていなければその目標地の蜜の量を-1します。

Step3：追従バチフェーズ

各追従バチが以下の確率に応じてルーレット選択で働きバチの目標地を決める。

\(p^i=f^i/\sum\nolimits_{n=1}^{N_{h}}f^{n}\)

\(i\)番目の働きバチの目標地に対するルーレット選択のための評価値\(f^i\)は以下のように計算します。

\(f^i=(vmax-v^i)/(vmax-vmin)\)

\(v^i\)は\(i\)番目の働きバチの目標地での目的関数の値です。目的関数の最小化では上の式を使います。なお\(vmax\)はすべての働きバチの目標地での目的関数の値のうちの最大値、\(vmin\)は最小値です。

ルーレット選択ではある乱数\(rand[0, 1]\)を発生させます。そして\(p^i\)を働きバチの順番で足していき、和が乱数を超えればそのときの働きバチの目的地を採用します。

つまり乱数が0.3のとき、1番目の働きバチの\(p^1\)が0.2、2番目の働きバチの\(p^2\)が0.15、なら2番目での和は0.35となり0.3を超えたので2番目の働きバチの目標地を採用します。

そして採用した目標値に対して働きバチフェーズ同様、新しい探索点\(\boldsymbol{y}^{i}\)を以下のように決定します。

\(y^i_j = \left\{ \begin{array}{l} x_j + rand[-1,1] \times ( x^i_j - x^l_j ) \quad j=k \\ x^i_j \qquad (otherwise) \end{array} \right. \)

改善されていなければその目標地の蜜の量を-1します。

Step4：偵察バチフェーズ

偵察バチが新しい目標地を見つけるフェーズです。

各働きバチの目標地での蜜の量が0になった目標地に対して、以下の式で新しい目標地を設定します。

(1)\(x_{j}^{i}=x_{j}^{L}+(x_{j}^{R}-x_{j}^{L})×rand[0, 1]\)

\(i\)番目の目標地のすべての成分に対して式(1)で新しい目標地を作成します。

目標地を作成したら蜜の量を\(h\)に戻します。

Step5：最良値の記録

現在のすべての働きバチの目標地での目的関数の値を計算し、その中の最良値が過去の最良値より優れていれば、より優れた目標地の座標に更新し、新たな最良値も更新します。

ABCアルゴリズムのプログラミング

ここからは実際にABCアルゴリズムの作成について書いていきます。

まずプログラムの全体をお示しして、その後それぞれの蜂のフェーズごとにプログラムの内容を解説していきます。

対象にしているのは次のようなSphere関数です。\(n\)は変数の数です。今回は\(n=3\)としました。

\(f(\boldsymbol{x})=\sum\nolimits_{i=1}^n x_i^2\)

ABCアルゴリズムのプログラム全体は以下となります。

import random
import copy
import math
xrange = []
for i in range(3):
    xrange.append([-5, 5])
harvest_bee = 7
follow_bee = 3
honey_max = 22
honey = []
flower = []
for i in range(harvest_bee):
    flower.append([])
    honey.append(honey_max)
    for j in range(len(xrange)):
        temp = xrange[j][0]+(xrange[j][1]-xrange[j][0])*random.random()
        if temp < xrange[j][0]:
            temp = xrange[j][0]
        if temp > xrange[j][1]:
            temp = xrange[j][1]
        flower[i].append(temp)        
iteration = 1000
#def calc(x):
#    sum1 = 0
#    sum2 = 0
#    for i in range(len(x)):
#        sum1 = sum1 + x[i]*x[i]
#        sum2 = sum2 + math.cos(2*math.pi*x[i])
#    return 20 - 20*math.exp(-0.2*math.sqrt((1/len(xrange))*sum1)) + math.exp(1) - math.exp((1/len(xrange))*sum2)
def calc(x):
    sum = 0
    for i in range(len(x)):
        sum = sum + x[i]*x[i]
    return sum
bestx = copy.copy(flower[0])
besty= calc(bestx)
def harvest():
    for i in range(harvest_bee):
        tempx = []
        l = random.randint(0, harvest_bee - 1)
        k = random.randint(0, len(xrange) - 1)
        for j in range(len(xrange)):
            if j == k:
                temp = flower[i][j] + random.uniform(-1, 1)*(flower[i][j]-flower[l][j])
                if temp < xrange[j][0]:
                    temp = xrange[j][0]
                if temp > xrange[j][1]:
                    temp = xrange[j][1]
                tempx.append(temp)                
            else:
                tempx.append(flower[i][j])
        tempy = calc(tempx)
        if tempy < calc(flower[i]):
            flower[i] = copy.copy(tempx)
            honey[i] = honey_max
        else:
            honey[i] = honey[i] - 1
def follow():
    sum = 0
    tempmin = 0
    tempmax = 0
    for i in range(harvest_bee):
        tempy = calc(flower[i])
        if tempy > calc(flower[tempmax]):
            tempmax = i
        if tempy < calc(flower[tempmin]):
            tempmin = i
    tempymax = calc(flower[tempmax])
    tempymin = calc(flower[tempmin])
    for i in range(harvest_bee):
        if tempymax != tempymin:
            f = (tempymax - calc(flower[i]))/(tempymax - tempymin)
            sum = sum + f
        else:
            f = 1
            sum = sum + f
    for i in range(follow_bee):
        p = random.random()
        sump = 0
        for j in range(harvest_bee):
            if tempymax != tempymin:
                f = (tempymax - calc(flower[j]))/(tempymax - tempymin)
                sump = sump + (f / sum)
            else:
                f = 1
                sump = sump + (f / sum)
            sump = sump + (f / sum)
            if sump > p:
                tempx = []
                l = random.randint(0, harvest_bee - 1)
                k = random.randint(0, len(xrange) - 1)
                for t in range(len(xrange)):
                    if t == k:
                        temp = flower[j][k] + random.uniform(-1, 1)*(flower[j][k]-flower[l][k])
                        if temp < xrange[t][0]:
                            temp = xrange[t][0]
                        if temp > xrange[t][1]:
                            temp = xrange[t][1]
                        tempx.append(temp)
                    else:
                        tempx.append(flower[j][t])
                tempy = calc(tempx)
                if tempy < calc(flower[j]):
                    flower[j] = copy.copy(tempx)
                    honey[j] = honey_max
                else:
                    honey[j] = honey[j] - 1
                break
def scout():
    for i in range(harvest_bee):
        if honey[i] <= 0:
            for j in range(len(xrange)):
                temp = xrange[j][0]+(xrange[j][1]-xrange[j][0])*random.random()
                if temp < xrange[j][0]:
                    temp = xrange[j][0]
                if temp > xrange[j][1]:
                    temp = xrange[j][1]
                flower[i][j] = temp
            honey[i] = honey_max
def record(bestx, besty):
    for i in range(harvest_bee):
        tempy = calc(flower[i])
        if tempy < besty:
            for j in range(len(xrange)):
                bestx[j] = flower[i][j]
    return calc(bestx)

for i in range(iteration):
    harvest()
    follow()
    scout()
    besty = record(bestx, besty)
print('best:', bestx)
print('besty:', besty)

Step1:初期化

初期化のプログラムは以下の部分です。

xrange = []
for i in range(3):
    xrange.append([-5, 5])
harvest_bee = 7
follow_bee = 3
honey_max = 22
honey = []
flower = []
for i in range(harvest_bee):
    flower.append([])
    honey.append(honey_max)
    for j in range(len(xrange)):
        temp = xrange[j][0]+(xrange[j][1]-xrange[j][0])*random.random()
        if temp < xrange[j][0]:
            temp = xrange[j][0]
        if temp > xrange[j][1]:
            temp = xrange[j][1]
        flower[i].append(temp)        
iteration = 1000

harvest_beeが働きバチの数、follow_beeが追従バチの数、honey_maxが蜜の上限です。

temp = xrange[j][0]+(xrange[j][1]-xrange[j][0])*random.random()

ここで新しい目標地を作成しています。

        if temp < xrange[j][0]:
            temp = xrange[j][0]
        if temp > xrange[j][1]:
            temp = xrange[j][1]

この部分で新しい目標地が変数の上下限を超えたらそれ以上いかないようにしています。これがないとどんどん変数の上下限を超える範囲を探索してしまい、収拾がつかなくなることがあります。

Step2:働きバチフェーズ

働きバチのフェーズは以下のようなプログラムになります。

def harvest():
    for i in range(harvest_bee):
        tempx = []
        l = random.randint(0, harvest_bee - 1)
        k = random.randint(0, len(xrange) - 1)
        for j in range(len(xrange)):
            if j == k:
                temp = flower[i][j] + random.uniform(-1, 1)*(flower[i][j]-flower[l][j])
                if temp < xrange[j][0]:
                    temp = xrange[j][0]
                if temp > xrange[j][1]:
                    temp = xrange[j][1]
                tempx.append(temp)                
            else:
                tempx.append(flower[i][j])
        tempy = calc(tempx)
        if tempy < calc(flower[i]):
            flower[i] = copy.copy(tempx)
            honey[i] = honey_max
        else:
            honey[i] = honey[i] - 1

以下の部分は働きバチが探索する部分です。

temp = flower[i][j] + random.uniform(-1, 1)*(flower[i][j]-flower[l][j])

良い目標地が新たに見つかった場合、見つからなかった場合それぞれの蜜の処理は以下の部分です。

        if tempy < calc(flower[i]):
            flower[i] = copy.copy(tempx)
            honey[i] = honey_max
        else:
            honey[i] = honey[i] - 1

Step3:追従バチフェーズ

追従バチの部分は以下となります。ルーレット選択の処理が入るので、結構長めとなっています。

def follow():
    sum = 0
    tempmin = 0
    tempmax = 0
    for i in range(harvest_bee):
        tempy = calc(flower[i])
        if tempy > calc(flower[tempmax]):
            tempmax = i
        if tempy < calc(flower[tempmin]):
            tempmin = i
    tempymax = calc(flower[tempmax])
    tempymin = calc(flower[tempmin])
    for i in range(harvest_bee):
        if tempymax != tempymin:
            f = (tempymax - calc(flower[i]))/(tempymax - tempymin)
            sum = sum + f
        else:
            f = 1
            sum = sum + f
    for i in range(follow_bee):
        p = random.random()
        sump = 0
        for j in range(harvest_bee):
            if tempymax != tempymin:
                f = (tempymax - calc(flower[j]))/(tempymax - tempymin)
                sump = sump + (f / sum)
            else:
                f = 1
                sump = sump + (f / sum)
            sump = sump + (f / sum)
            if sump > p:
                tempx = []
                l = random.randint(0, harvest_bee - 1)
                k = random.randint(0, len(xrange) - 1)
                for t in range(len(xrange)):
                    if t == k:
                        temp = flower[j][k] + random.uniform(-1, 1)*(flower[j][k]-flower[l][k])
                        if temp < xrange[t][0]:
                            temp = xrange[t][0]
                        if temp > xrange[t][1]:
                            temp = xrange[t][1]
                        tempx.append(temp)
                    else:
                        tempx.append(flower[j][t])
                tempy = calc(tempx)
                if tempy < calc(flower[j]):
                    flower[j] = copy.copy(tempx)
                    honey[j] = honey_max
                else:
                    honey[j] = honey[j] - 1
                break

まずルーレット選択するために最良値と最悪値を求めます。

    tempmin = 0
    tempmax = 0
    for i in range(harvest_bee):
        tempy = calc(flower[i])
        if tempy > calc(flower[tempmax]):
            tempmax = i
        if tempy < calc(flower[tempmin]):
            tempmin = i
    tempymax = calc(flower[tempmax])
    tempymin = calc(flower[tempmin])

次にルーレット選択するための総和を求めます。

    for i in range(harvest_bee):
        if tempymax != tempymin:
            f = (tempymax - calc(flower[i]))/(tempymax - tempymin)
            sum = sum + f
        else:
            f = 1
            sum = sum + f

なぜtempymax != tempymin:の条件分岐があるのか、という点ですが、働きバチが少なかったり、値が収束してくると、目的関数がすべて横並びの状態になって、tempymax – tempyminがゼロになることがあるので、ゼロで割らないために横並びならすべての働きバチが選ばれる確率を同等としています。

次にルーレット選択で目的関数の値がよい目標地を優先的に選ぶ処理です。

    for i in range(follow_bee):
        p = random.random()
        sump = 0
        for j in range(harvest_bee):
            if tempymax != tempymin:
                f = (tempymax - calc(flower[j]))/(tempymax - tempymin)
                sump = sump + (f / sum)
            else:
                f = 1
                sump = sump + (f / sum)
            sump = sump + (f / sum)
            if sump > p:
                tempx = []
                l = random.randint(0, harvest_bee - 1)
                k = random.randint(0, len(xrange) - 1)
                for t in range(len(xrange)):
                    if t == k:
                        temp = flower[j][k] + random.uniform(-1, 1)*(flower[j][k]-flower[l][k])
                        if temp < xrange[t][0]:
                            temp = xrange[t][0]
                        if temp > xrange[t][1]:
                            temp = xrange[t][1]
                        tempx.append(temp)
                    else:
                        tempx.append(flower[j][t])
                tempy = calc(tempx)
                if tempy < calc(flower[j]):
                    flower[j] = copy.copy(tempx)
                    honey[j] = honey_max
                else:
                    honey[j] = honey[j] - 1
                break

sumpはどんどん各働きバチの目的関数の値によって増えていきます。目的関数が良ければ大きく増えるし、良くなければほとんど増えません。sumpが乱数pを超えたらそのときの働きバチを選択します。

働きバチを選択したら、あとは働きバチフェーズと同様探索を行います。

Step4:偵察バチフェーズ

各働きバチの目標地の蜜の量に応じて新しい目標地を作成します。

    for i in range(harvest_bee):
        if honey[i] <= 0:
            for j in range(len(xrange)):
                temp = xrange[j][0]+(xrange[j][1]-xrange[j][0])*random.random()
                if temp < xrange[j][0]:
                    temp = xrange[j][0]
                if temp > xrange[j][1]:
                    temp = xrange[j][1]
                flower[i][j] = temp
            honey[i] = honey_max

Step5:最良値の記録

現在の最良値を保存します。

    for i in range(harvest_bee):
        tempy = calc(flower[i])
        if tempy < besty:
            for j in range(len(xrange)):
                bestx[j] = flower[i][j]
    return calc(bestx)

return calc(bestx)の部分は関数の引数にbestyを設定しても参照渡しできないので、いったん出力して元の値に代入することで値を更新します。

besty = record(bestx, besty)

この部分です。

ABCアルゴリズムの優れた点【遺伝的アルゴリズムとの比較】

遺伝的アルゴリズムはランダムサーチの突然変異でしか細かい値を見つけられない

遺伝的アルゴリズムでは交叉で良い値を選別していくため、連続変数になるとより細かい値を重点的に探す機能がいまいちです。細かい新しい探索点は突然変異のランダムサーチでしか見つからないので、なかなか狙った位置に突然変異した値が来ません。摂動という突然変異である程度解消できますが、摂動幅をどうするのかというのは悩ましい問題です。

ABCアルゴリズムは以下の式で重点的な探索を可能としています。

\(y^i_j = \left\{ \begin{array}{l} x_j + rand[-1,1] \times ( x^i_j - x^l_j ) \quad j=k \\ x^i_j \qquad (otherwise) \end{array} \right. \)

新しい目標地が生成される最大幅がすべての目標地の一番遠いところ同士から選択されるので、全体が集まってくると最大幅も小さくなります。これで重点的な探索が可能です。

ABCアルゴリズムは局所解に陥りにくい

ABCアルゴリズムは偵察バチフェーズで蜜を取りつくした見込みのない目標地をリセットする機能があります。

蜜を取りつくすということは、その目標地近辺で探索しても良い値が見つからなかったということです。そういうときに目標地をリセットしてランダムサーチすることで局所解以外の場所を見つけやすくしています。

これは遺伝的アルゴリズムの突然変異に相当しますが、遺伝的アルゴリズムが本当にランダムに突然変異させていたところ、ABCアルゴリズムは見込みのない目標地の働きバチだけ突然変異させているため、見込みのある目標地の探索が突然変異で打ち切られにくいというところが優れていますね。

遺伝的アルゴリズムおよびソルバーのエボリューショナリーとの比較

せっかくなので、過去に当ブログで進化計算で計算した結果と同じような条件でベンチマーク問題を解いてみました。

まずは以下のブログのSphere関数です。この記事では3変数です。

遺伝的アルゴリズムでの結果は繰り返し1000回で探索点10個の場合以下となりました。

[0.016800221679470262, 0.014557315436035445, 0.0013207670648451852]

0.0004959073068231582

同じ条件に近づけて以下のパラメータでABCアルゴリズムを実行した結果は以下となりました。

働きバチ：6
追従バチ：4
蜜の最大量：20
繰り返し回数：1000

ABCアルゴリズムの計算結果は以下となりました。

best: [-2.9133385789604886e-17, 1.7713002753217187e-17, -3.1047929529006967e-15]

besty: 9.64090178501593e-30

続いては以下の記事でのエクセルソルバーのエボリューショナリーとの結果比較です。

対象としたのは2変数のAckley関数です。

エボリューショナリーの結果は

[4.33e-06, 7.58e-07]

1.24e-05

でした。ABCアルゴリズムの結果はパラメータを上のSphere関数と同じにして以下となりました。

best: [-3.261885767616958e-13, 1.8510243215913983e-13]

besty: 1.0622613899613498e-12

どちらのベンチマーク問題でも過去の計算結果より良い値が求まっています。ABCアルゴリズムは結構優れた手法のようですね。

パラメータの調整

ABCアルゴリズムではパラメータをどうするかという問題があります。特に働きバチと追従バチの数のバランスと、蜜の最大量は探索性能を大きく変動させるパラメータです。

最適なパラメータを決めるのは一種の最適化問題となっています。

そこでせっかくABCアルゴリズムを作成したのだから、ABCアルゴリズムでABCアルゴリズムの最適なパラメータを探索してみようと思い立ちました。

ここからはご興味のある方向けとなります。プログラムには継承などを利用しますが、プログラミングの手法については詳しく解説しません。ある程度プログラミングができることを前提に書いていきます。一からすべて解説すると非常に長くなってしまうからです。

では早速そのプログラムをお示しします。

import random
import copy
import math
class abc:
    def __init__(self, u_l, harvest_bee, follow_bee, honey_max, iteration):
        self.xrange = u_l
        self.harvest_bee = int(math.ceil(harvest_bee))
        self.follow_bee = int(math.ceil(follow_bee))
        self.honey_max = int(math.ceil(honey_max))
        self.honey = []
        self.flower = []
        self.flower_y = []
        for i in range(self.harvest_bee):
            self.flower.append([])
            self.honey.append(self.honey_max)
            for j in range(len(self.xrange)):
                temp = self.xrange[j][0]+(self.xrange[j][1]-self.xrange[j][0])*random.random()
                if temp < self.xrange[j][0]:
                    temp = self.xrange[j][0]
                if temp > self.xrange[j][1]:
                    temp = self.xrange[j][1]
                self.flower[i].append(temp)
            self.flower_y.append(self.calc(self.flower[i]))
        self.iteration = iteration
        self.bestx = copy.copy(self.flower[0])
        self.besty= self.flower_y[0]
    def calc(self, x):
        sum1 = 0
        sum2 = 0
        for i in range(len(x)):
            sum1 = sum1 + x[i]*x[i]
            sum2 = sum2 + math.cos(2*math.pi*x[i])
        return 20 - 20*math.exp(-0.2*math.sqrt((1/len(self.xrange))*sum1)) + math.exp(1) - math.exp((1/len(self.xrange))*sum2)
#    def calc(self, x):
#        sum = 0
#        for i in range(len(x)):
#            sum = sum + x[i]*x[i]
#        return sum
    def harvest(self):
        for i in range(self.harvest_bee):
            tempx = []
            l = random.randint(0, self.harvest_bee - 1)
            k = random.randint(0, len(self.xrange) - 1)
            for j in range(len(self.xrange)):
                if j == k:
                    temp = self.flower[i][j] + random.uniform(-1, 1)*(self.flower[i][j]-self.flower[l][j])
                    if temp < self.xrange[j][0]:
                        temp = self.xrange[j][0]
                    if temp > self.xrange[j][1]:
                        temp = self.xrange[j][1]
                    tempx.append(temp)
                else:
                    tempx.append(self.flower[i][j])
            tempy = self.calc(tempx)
            if tempy < self.flower_y[i]:
                self.flower[i] = copy.copy(tempx)
                self.honey[i] = self.honey_max
                self.flower_y[i] = tempy
            else:
                self.honey[i] = self.honey[i] - 1
    def follow(self):
        sum = 0
        tempmin = 0
        tempmax = 0
        for i in range(self.harvest_bee):
            tempy = self.flower_y[i]
            if tempy > self.flower_y[tempmax]:
                tempmax = i
            if tempy < self.flower_y[tempmin]:
                tempmin = i
        tempymax = self.flower_y[tempmax]
        tempymin = self.flower_y[tempmin]
        for i in range(self.harvest_bee):
            if tempymax - tempymin != 0.0:
                f = (tempymax - self.flower_y[i])/(tempymax - tempymin)
                sum = sum + f
            else:
                f = 1
                sum = sum + f
        for i in range(self.follow_bee):
            p = random.random()
            sump = 0
            for j in range(self.harvest_bee):
                if tempymax - tempymin != 0.0:
                    f = (tempymax - self.flower_y[j])/(tempymax - tempymin)
                    sump = sump + (f / sum)
                else:
                    f = 1
                    sump = sump + (f / sum)
                if sump > p:
                    tempx = []
                    l = random.randint(0, self.harvest_bee - 1)
                    k = random.randint(0, len(self.xrange) - 1)
                    for t in range(len(self.xrange)):
                        if t == k:
                            temp = self.flower[j][k] + random.uniform(-1, 1)*(self.flower[j][k]-self.flower[l][k])
                            if temp < self.xrange[t][0]:
                                temp = self.xrange[t][0]
                            if temp > self.xrange[t][1]:
                                temp = self.xrange[t][1]
                            tempx.append(temp)
                        else:
                            tempx.append(self.flower[j][t])
                    tempy = self.calc(tempx)
                    if tempy < self.flower_y[j]:
                        self.flower[j] = copy.copy(tempx)
                        self.honey[j] = self.honey_max
                        self.flower_y[j] = tempy
                    else:
                        self.honey[j] = self.honey[j] - 1
                    break
    def scout(self):
        for i in range(self.harvest_bee):
            if self.honey[i] < 0:
                for j in range(len(self.xrange)):
                    temp = self.xrange[j][0]+(self.xrange[j][1]-self.xrange[j][0])*random.random()
                    if temp < self.xrange[j][0]:
                        temp = self.xrange[j][0]
                    if temp > self.xrange[j][1]:
                        temp = self.xrange[j][1]
                    self.flower[i][j] = temp
                self.honey[i] = self.honey_max
    def record(self, bestx, besty):
        for i in range(self.harvest_bee):
            tempy = self.flower_y[i]
            if tempy < besty:
                for j in range(len(self.xrange)):
                    bestx[j] = self.flower[i][j]
                besty = tempy
        return besty
    def run(self):
        for i in range(self.iteration):
            self.harvest()
            self.follow()
            self.scout()
            self.besty = self.record(self.bestx, self.besty)
class parametersearch(abc):
    def calc(self, x):
        xtemp = []
        xtemp.append(x[0])
        xtemp.append(10-x[0])
        xtemp.append(x[1])
        parent = abc([[-10, 10],[-10, 10]], xtemp[0], xtemp[1], xtemp[2], 500)
        parent.run()
        return parent.besty
abcrunsearch = parametersearch([[3, 8], [3, 30]], 5, 3, 10, 50)
abcrunsearch.run()
print('best:', abcrunsearch.bestx)
print('besty:', abcrunsearch.besty)

簡単に説明するとABCアルゴリズムのクラスを作成して、パラメータ最適化用のABCアルゴリズムを継承によって作成しています。関数の計算部分が呼び出しごとに一意に決まらないので一時的に保存しておくflower_yリストを導入しています。

継承したクラスの中で関数の計算部分だけ変更しています。計算部分で通常のABCアルゴリズムを実行して、その最良値を継承したクラスの目的関数の値にしています。

知りたいのは働きバチと追従バチの比なので、この二種類の蜂の数の合計は10としました。蜜の上限は適当に決めました。

実際はABCアルゴリズムが確率的な手法なので、ある変数に対して必ず一意に同じ値が出るわけではないので、出てきた最適なパラメータの信頼性は確実ではありません。

しかしより良い値を探すという手法で吐き出された最適なパラメータに何らかの傾向があれば、それはそれで有益な情報なのではないか、という方針のもと、初心者が実験的にやってみた程度の感覚で見ていただければと思います。

また初心者の実験なので、対象とする最適化問題も2変数のAckley関数と3変数のSphere関数のみです。

大学の研究レベルの信頼性のある実験ではありません。結果に対していかなる保証もしないという前提のもとご覧ください。

なおAckley functionは以下となります。

\({f(x_{1} \cdots x_{n})=20-20\exp \biggl( -0.2\sqrt{\frac{1}{n}\sum_{i=1}^{n}x_{i}^2} \biggr) +e-\exp \biggl(\frac{1}{n}\sum_{i=1}^{n}\cos(2\pi x_{i}) \biggr)
}\)

\({-10 \leq x_{i} \leq 10, f_{min}(0, \cdots , 0)=0}\)

パラメータチューニングの結果

3変数のSphere関数と2変数のAckley関数をそれぞれ10回パラメータチューニングしてみた結果は以下となりました。

まずはSphere関数。

働きバチの数	蜜の最大数
3.9764	28.72408
4.663546	24.22744
4.704776	29.19044
4.799971	27.44192
4.748477	30
4.074736	24.33931
4.696321	23.94735
4.164656	26.15571
4.384457	27.1224
4.246951	28.9874

平均値は

働きバチの数：4.446029

蜜の最大数：27.01361

続いてはAckley関数です。

働きバチの数	蜜の最大数
7.334909	24.79809
6.355864	19.41943
4.379705	24.9187
6.133522	23.6597
6.320366	20.37498
8	29.59223
7.184658	17.27927
7.327894	21.17496
7.046137	16.54497
7.987263	26.36477