Optunaで多目的最適化を実行

import optuna
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from optuna.pruners import SuccessiveHalvingPruner

# 目的関数を定義する関数
def VLMOP2(x):
    n = len(x)
    f1 = -(1 - np.exp(-1 * np.sum((x - 1/np.sqrt(n)) ** 2)))
    f2 = -(1 - np.exp(-1 * np.sum((x + 1/np.sqrt(n)) ** 2)))
    return f1, f2

#探索点を評価する関数
def objective(trial):
    n = 2  # 関数の次元数
    x = [trial.suggest_float(f'x{i}', -5, 5, step=0.01) for i in range(n)]  # 探索範囲を-5から5まで1000刻みに設定
    return VLMOP2(x)

#VLMOP2のPareto frontの計算
x = np.linspace(-5, 5, 100)
X1, X2 = np.meshgrid(x, x)
points = np.column_stack((X1.ravel(), X2.ravel()))
y = np.array([VLMOP2(point) for point in points])
pareto_front_vlmop2 = []
for point in y:
    is_pareto = True
    for other_point in y:
        if (other_point[0] >= point[0] and other_point[1] >= point[1] and
            (other_point[0] > point[0] or other_point[1] > point[1])):
            is_pareto = False
            break
    if is_pareto:
        pareto_front_vlmop2.append(point)
pareto_front_vlmop2 = np.array(pareto_front_vlmop2)

# 多目的最適化及び探索進捗を可視化する関数
def optimize_and_plot(num_trials, num_frames):
    study = optuna.create_study(directions=['maximize', 'maximize'], sampler=optuna.samplers.MOTPESampler(),
                                pruner=SuccessiveHalvingPruner())
    animation_data = []

    def update(frame):
        global pareto_front
        study.optimize(objective, n_trials=1)  # フレームごとのトライアル数を1に変更
        trials = study.trials_dataframe()
        f1, f2 = trials['values_0'], trials['values_1']
        ax.clear()
        ax.scatter(f1, f2, label='Optimized points')  # 凡例を追加
        ax.set_xlabel('f1')
        ax.set_ylabel('f2')
        ax.set_title(f'Optimization Progress (Frame {frame+1}/{num_frames})')
        
        # Pareto frontの解を赤色でハイライト表示
        pareto_front = np.array([t.values for t in study.best_trials])  # 最良トライアルからPareto frontを取得
        ax.scatter(pareto_front[:, 0], pareto_front[:, 1], color='red', marker='o', label='Pareto front')
        
        # VLMOP2のPareto frontを黒色でプロット
        ax.plot(pareto_front_vlmop2[:, 0], pareto_front_vlmop2[:, 1], color='black', linewidth=2, label='VLMOP2 Pareto front')

        ax.legend()  # 凡例を表示

    fig, ax = plt.subplots()
    ani = FuncAnimation(fig, update, frames=num_frames, repeat=False)
    ani.save('multi_objective_optimization.gif', writer='pillow', fps=7)

#パレート超体積を計算する関数
def pareto_hypervolume(pareto_front):
    
    # 探索点のパレートフロントのf1を昇順に並び替え
    pareto_front = pareto_front[np.argsort(pareto_front[:, 0])]
    
    # パレート超体積の初期化
    hypervolume = 0.0
    
    # パレート超体積の計算
    for i in range(len(pareto_front) - 1):
        width = pareto_front[i + 1, 0] - pareto_front[i, 0]
        height = pareto_front[i, 1]
        hypervolume += width * (height+1)
    
    return hypervolume    
    
#多目的最適化を実行
optimize_and_plot(num_trials=200, num_frames=200)

# パレート超体積の計算を実行
hypervolume = pareto_hypervolume(pareto_front)
vlmop2_hypervolume = pareto_hypervolume(pareto_front_vlmop2)
print(f'VLMOP2のパレートフロント : {vlmop2_hypervolume}')
print(f'探索点のパレートフロント : {hypervolume}')

'''
VLMOP2のパレートフロント : 0.34549703030586165
探索点のパレートフロント : 0.31021040536297584
'''

def VLMOP2(x):
    n = len(x)
    f1 = -(1 - np.exp(-1 * np.sum((x - 1/np.sqrt(n)) ** 2)))
    f2 = -(1 - np.exp(-1 * np.sum((x + 1/np.sqrt(n)) ** 2)))
    return f1, f2

#探索点を評価する関数
def objective(trial):
    n = 2  # 関数の次元数
    x = [trial.suggest_float(f'x{i}', -5, 5, step=0.01) for i in range(n)]  # 探索範囲を-5から5まで1000刻みに設定
    return VLMOP2(x)

#VLMOP2のPareto frontの計算
x = np.linspace(-5, 5, 100)
X1, X2 = np.meshgrid(x, x)
points = np.column_stack((X1.ravel(), X2.ravel()))
y = np.array([VLMOP2(point) for point in points])
pareto_front_vlmop2 = []
for point in y:
    is_pareto = True
    for other_point in y:
        if (other_point[0] >= point[0] and other_point[1] >= point[1] and
            (other_point[0] > point[0] or other_point[1] > point[1])):
            is_pareto = False
            break
    if is_pareto:
        pareto_front_vlmop2.append(point)
pareto_front_vlmop2 = np.array(pareto_front_vlmop2)

# 多目的最適化及び探索進捗を可視化する関数
def optimize_and_plot(num_trials, num_frames):
    study = optuna.create_study(directions=['maximize', 'maximize'], sampler=optuna.samplers.MOTPESampler(),
                                pruner=SuccessiveHalvingPruner())
    animation_data = []

    def update(frame):
        global pareto_front
        study.optimize(objective, n_trials=1)  # フレームごとのトライアル数を1に変更
        trials = study.trials_dataframe()
        f1, f2 = trials['values_0'], trials['values_1']
        ax.clear()
        ax.scatter(f1, f2, label='Optimized points')  # 凡例を追加
        ax.set_xlabel('f1')
        ax.set_ylabel('f2')
        ax.set_title(f'Optimization Progress (Frame {frame+1}/{num_frames})')
        
        # Pareto frontの解を赤色でハイライト表示
        pareto_front = np.array([t.values for t in study.best_trials])  # 最良トライアルからPareto frontを取得
        ax.scatter(pareto_front[:, 0], pareto_front[:, 1], color='red', marker='o', label='Pareto front')
        
        # VLMOP2のPareto frontを黒色でプロット
        ax.plot(pareto_front_vlmop2[:, 0], pareto_front_vlmop2[:, 1], color='black', linewidth=2, label='VLMOP2 Pareto front')

        ax.legend()  # 凡例を表示

    fig, ax = plt.subplots()
    ani = FuncAnimation(fig, update, frames=num_frames, repeat=False)
    ani.save('multi_objective_optimization.gif', writer='pillow', fps=7)

#パレート超体積を計算する関数
def pareto_hypervolume(pareto_front):
    
    # 探索点のパレートフロントのf1を昇順に並び替え
    pareto_front = pareto_front[np.argsort(pareto_front[:, 0])]
    
    # パレート超体積の初期化
    hypervolume = 0.0
    
    # パレート超体積の計算
    for i in range(len(pareto_front) - 1):
        width = pareto_front[i + 1, 0] - pareto_front[i, 0]
        height = pareto_front[i, 1]
        hypervolume += width * (height+1)
    
    return hypervolume    

#多目的最適化を実行
optimize_and_plot(num_trials=200, num_frames=200)

# パレート超体積の計算を実行
hypervolume = pareto_hypervolume(pareto_front)
vlmop2_hypervolume = pareto_hypervolume(pareto_front_vlmop2)
print(f'VLMOP2のパレートフロント : {vlmop2_hypervolume}')
print(f'探索点のパレートフロント : {hypervolume}')