🐙

Optunaを使ってランベール問題を最適化してみた

2023/04/03に公開

Optuna

軌道力学

tech

はじめに

はじめまして！初投稿になります．普段は某企業でエンジニアしており，宇宙系プロジェクトに携わらせてもらってます．リスキリングの一環でインプット/アウトプットをちょくちょくしていこうかなと思い始めました．
今回は機械学習のハイパーパラメータ調整などに利用されるOptunaを使ってみたいと思う．
よく最適化関連の例にはHimmelblau関数がサンプルとして用いられるが，宇宙系エンジニアとしては軌道力学におけるランベール問題を題材にして，Optunaを使う練習をしてみようと思う（Optunaを使うまでもないが）．

Optunaについて

Optunaは，ブラックボックス最適化フレームワークのひとつである．グリッドサーチやランダムサーチだけでなく，ベイズ最適化や進化計算，準モンテカルロ法などのアルゴリズムを提供しており，非常に有用なツールになっている．さらには大規模な機械学習モデルを学習する場合などに，RDB(Relational Database)を利用した分散並列最適化もできる．

ランベール問題について

地球から火星へ遷移する際の軌道を題材にして，この軌道をランベール問題^[1]を用いて簡易的に解いてみよう．ランベール問題は，2点（出発地点と到着地点）と遷移時間が与えられたとき，出発時の速度と到着時の速度が求まるというもの．このランベール問題を利用して地球出発時の脱出相対速度ベクトル（色付き等高線，図ではC3になっている）と火星到着時の接近相対速度ベクトル（青等高線）（が最小となる出発/到着日）を最適化してみようと思う．以下の図はポークチョップ図と呼ばれている．

使用したコードは以下のリポジトリにあります．

ここで地球と火星の位置はSPICE/spiceypyを，ランベール問題のソルバにはpoliastroを利用した．

単目的最適化

まずはシンプルに地球出発時の脱出相対速度のみを最適化する．Optunaでは1回の評価をトライアルと呼び，トライアルを繰り返す一連の最適化プロセスをスタディと呼ぶ．目的関数の引数にTrialオブジェクトを渡す．suggest APIというメソッドを利用して目的関数への入力変数を提案する．

def objective(trial):
    # suggest API
    x = trial.suggest_float("x", dep_time1, dep_time2)
    y = trial.suggest_float("y", arr_time1, arr_time2)

ここで，dep_time1, dep_time2, arr_time1, arr_time2はそれぞれ出発/到着日時の探索範囲を設定しており，今回はそれぞれをMJD形式の実数型としてsuggest APIに入力した．objective関数内でランベール問題を解き，地球出発時の脱出相対速度を返す．

Studyオブジェクトを生成する．Studyオブジェクト生成時に，最小化するか最大化するか指定，サンプラーを指定する．デフォルトではサンプラーはTPESampler(tree-structured Parzen estimator)が選ばれる．optimizeメソッドに上記で定義したobjective関数，ループ回数を入力し最適化を実行する．

    # default, TPE(ベイズ最適化)
    study = optuna.create_study(direction="minimize")
    # optimize
    study.optimize(objective, n_trials=100)

    print(f"Best objective value: {study.best_value}")
    print(f"Best parameter: {study.best_params}")
    dep_opt = astropy.time.Time(study.best_params["x"], format='mjd')
    arr_opt = astropy.time.Time(study.best_params["y"], format='mjd')
    print(dep_opt.datetime)
    print(arr_opt.datetime)

以下の実行結果が得られる．

Best objective value: 11.141860022519491
Best parameter: {'x': 60587.7370035073, 'y': 60936.865608620465}
2024-10-04 17:41:17.103031
2025-09-18 20:46:28.584808

ポークチョップ図上で見ると，たしかに最適解が求まっていることがわかる．

デフォルトではサンプラーはTPEだが，以下のようにサンプラーを明示的にグリッドサーチなどに変更することもできる．サンプラーにはTPESampler（ベイズ最適化），NSGAIISampler（進化計算），BoTorchSampler（ベイズ最適化），QMCSampler（準モンテカルロ法），RandomSampler（ランダムサーチ），GridSampler（グリッドサーチ）などが用意されている．

    # グリッドサーチ
    search_space = {
        'x': np.linspace(dep_time1, dep_time2, 10), 
        'y': np.linspace(arr_time1, arr_time2, 10)
    }
    study = optuna.create_study(
        direction="minimize", 
        sampler=optuna.samplers.GridSampler(search_space=search_space)
    )

また，DBを使用して学習を保存することもできる．

    #
    # RDB
    #
    study = optuna.create_study(
        direction="minimize", 
        storage="sqlite:///optuna.db", 
        study_name="ch1-conditional", 
    )
    #
    # load study
    #
    study = optuna.load_study(
        storage="sqlite:///optuna.db", 
        study_name="ch1-conditional", 
    )

多目的最適化

次に地球出発時の脱出相対速度と火星到着時の接近相対速度を同時に最適化する．変更点はobjective関数内で，地球出発時の脱出相対速度と火星到着時の接近相対速度を返すことと，Studyオブジェクト生成時に目的ごとの最適化方向を指定することである．
最適化結果を表示するコードを以下のように変更した．単目的の時とは異なり，最良のトライアルが一意には定まらず，複数のトライアルが表示される．

    print("Best Trials")
    for t in study.best_trials:
        print(f"- [{t.number}] params={t.params}, values={t.values}")
        dep_opt = astropy.time.Time(t.params["x"], format='mjd')
        arr_opt = astropy.time.Time(t.params["y"], format='mjd')
        print(dep_opt.datetime)
        print(arr_opt.datetime)

多目的最適化の場合には，最良のトライアルが一意には定まらず，トレードオフの関係にあることが多い．トライアル群をパレートフロントと呼ばれ，Optunaではパレートフロントを可視化することができる．

    fig = optuna.visualization.plot_pareto_front(
        study, 
        include_dominated_trials=False
    )
    fig.write_image('./img/pareto_front.png')  # save as png file ($ pip install -U kaleido)

798回目のトライアルでは地球出発時の脱出相対速度が最小となり，958回目のトライアルでは火星到着時の接近相対速度が最小となる．

この結果をポークチョップ図に合わせてみると，多目的最適化が行えていることがわかる．

制約つき最適化

今回のように目的が複数（地球出発時の脱出相対速度と火星到着時の接近相対速度の最小化）ある場合，むやみに多目的最適化を採用するよりも，最適化対象の目的を1つに絞り，それ以外を制約として表現することができる場合には制約あり最適化を採用する方が良いこともある．
ここでは，火星到着時の接近相対速度を2.46[km/s]以下に抑えつつ，地球出発時の脱出相対速度を最小化してみる．
objective関数の変更点は「火星到着時の接近相対速度が2.46以下」という制約をTrialオブジェクトのユーザ属性として設定している点である．

    # constraints
    trial.set_user_attr("constraints", [RANGE_arr - 2.46])

またサンプラーオブジェクト生成時に，制約関数を設定する．

    # サンプラーに制約をもたせる
    sampler = optuna.samplers.TPESampler(
        constraints_func=lambda trial: trial.user_attrs["constraints"]
    )

以下のような実行結果が得られた．

- [470] params={'x': 60586.76392099191, 'y': 60917.95700574684}, values=[3.3570775049023354, 2.4585219227303168]
2024-10-03 18:20:02.773701
2025-08-30 22:58:05.296527

この解はポークチョップ図で書くと以下のようなイメージである．

おわりに

とまぁこのような感じでOptunaに触れてみることができた．今後なにかよい使い道を思いついたらアウトプットしていきたい．

参考文献

Optunaによるブラックボックス最適化, オーム社

脚注

"Revisiting Lambert's problem", Dario Izzo ↩︎