前回まではgradientベースの最適化を行った。解析的な勾配であれ、数値的な（差分法）勾配であれ、目的関数の変化量から最適な方向を探索することになる。この方法は、計算時間は短くてすむが、目的関数の変化が複雑（多峰性を持つ）な場合、必ずしも大域的な（ベストな）最適値に収束するとは限らず、初期値に近い局所的な（ベターな）最適値に落ち着いてしまう可能性がある。従って、勾配ベースの手法は、現状からの「改良」には向いているが、様々なな設計変数の範囲から、新たに「ベスト」な設計変数の組を見出すには不向きと考えられる。
このような場合を防ぐ大域的な方法として、勾配に依存しない手法（derivertive-free method）が提案されている。古くは、モンテカルロ法のような「しらみつぶし型」があるが、近年は生物の進化過程をモデル化した手法がしばしば用いられ、EA（Evolutional Algorithm）として総称されているようである。
Dakotaで導入されているパッケージ名としては(括弧内は手法名)、COLINY(PS, EA,...), JEGA(SOGA/MOGA),EGO(?),DIIRECT(?),OPT++(?)などがある。今回は、JEGAパッケージを使った多目的関数の最適化を考えてみる。
minimize: f1=2*sqrt(x1), f2=x1-x1*x2+5
subject to 1<=x1<=4, 1<=x2<=2
<dakota.in>
変更部分はmethodセクションのみである。大分、記述が増えている。
environment
graphics
tabular_graphics_data
tabular_graphics_file = 'objective.dat'
method
moga // MOGA(多目的GA)
output silent //画面アウトプットコントロール
seed = 10983 // 乱数発生用シード
final_solutions = 3 // 最終的なBest solutionとして残す設計変数組の数
max_function_evaluations = 2500 // 目的関数の最大評価回数
initialization_type unique_random // 初期値生成方法
crossover_type shuffle_random // 交配(crossover)の種類
num_offspring = 2 num_parents = 2 //
crossover_rate = 0.8 //交配の確率
mutation_type replace_uniform //変異(mutation)の種類
mutation_rate = 0.1 //変異の確率
fitness_type domination_count //
replacement_type below_limit = 6 // 設計組をまとめる下限値
shrinkage_percentage = 0.9 // ?
convergence_type metric_tracker // ?
percent_change = 0.05 num_generations = 40 // 40世代
model
single
variables
continuous_design = 2
initial_point 2.0 1.5
lower_bounds 1.0 1.0
upper_bounds 4.0 2.0
descriptors 'x1' "x2"
interface
fork
analysis_driver = 'simulator_script'
parameters_file = 'params.in'
results_file = 'results.out'
work_directory directory_tag
copy_files = 'templatedir/*'
named 'workdir' file_save directory_save
aprepro
responses
objective_functions = 2
no_gradients
no_hessians
以上の設定で実施してみると、以下の通りとなる。パレートフロントが明確に出ているが、最大評価回数2500回をすべて計算しているため、計算時間が非常にかかる。とくに、ソルバーがCFDやFEAの計算の場合には、現実的な計算時間とは言えなくなる。そこで、いくつか工夫を凝らす必要がある。
①並列計算の対応
設計変数の組の1個1個は独立に評価されるので、マルチコアの計算機であれば、以下のinterfaceの部分に2行追加することで、非同期(asyncribous)で、同時に目的関数の評価を４つ行うことになる。
interface
[snip]
asynchronousevaluation_concurrency = 4
[snip]
※ プロセスを見ると、４つ起動していることが確認できるが、いろいろオーバーヘッドもあり、必ずしも計算時間が1/4にはならないようです。
この他、dakotaそのもののプロセスを並列実行もできるようである（未確認）。
mpirun -np 2 dakota -i dakota.in
②Surrogateモデルによる応答曲面の近似（surrogate based optimaization）
設計変数の組を根本的に減らすため、得られた目的関数の応答曲面を近似（surrogated model）し、新たな設計変数の組を生成し、目的関数を評価する。
SBOの基本的な流れは、下記の通りである。
- DACEによる設計変数のdata sampling
- サンプリングした設計変数による目的関数の評価（e.g., CFD/FEA simulation）
- メタモデリング（目的関数評価値を使った近似(surrogated)応答曲面fittingモデルの作成）
- 最適化の実施
このメタモデリングの際、幾つかの設計変数の組では実際の目的関数が得られていることから、これらのデータを補間して近似面を形成する。補間方法には、最小自乗法やクリギング法、多項式近似などがある。
下記ように書き換えたサンプルを、kriging+GAで実行してみた。
※ -oオプションで出力を指定していないことに注意。指定するとなぜか実行されない。
<dakota_sbo.in>
# Dakota Input File: dakota_sbo.in
# Usage:
# dakota -i dakota_sbo.in > dakota.stdout
environment
# graphics
tabular_graphics_data
tabular_graphics_file = 'objective.dat'
method_pointer = 'SBO'
method
id_method = 'SBO'
surrogate_based_global
model_pointer = 'SURROGATE'
approx_method_pointer = 'MOGA'
max_iterations = 5
replace_points
method
id_method = 'MOGA'
moga
output silent
seed = 10983
population_size = 300
max_function_evaluations = 2500
initialization_type unique_random
crossover_type shuffle_random
num_offspring = 2 num_parents = 2
crossover_rate = 0.8
mutation_type replace_uniform
mutation_rate = 0.1
fitness_type domination_count
replacement_type below_limit = 6
shrinkage_percentage = 0.9
niching_type distance 0.05 0.05
postprocessor_type
orthogonal_distance 0.05 0.05
convergence_type metric_tracker
percent_change = 0.05 num_generations = 10
model
id_model = 'SURROGATE'
surrogate global
dace_method_pointer = 'DACE'
correction additive first_order
gaussian_process dakota
method
id_method = 'DACE'
model_pointer = 'TRUTH'
sampling
samples = 30 seed = 531
sample_type lhs
model
id_model = 'TRUTH'
single
variables
continuous_design = 2
initial_point 2.0 1.5
lower_bounds 1.0 1.0
upper_bounds 4.0 2.0
descriptors 'x1' "x2"
interface
fork
# asynchronous
# evaluation_concurrency = 4
analysis_driver = 'simulator_script'
parameters_file = 'params.in'
results_file = 'results.out'
work_directory directory_tag
copy_files = 'templatedir/*'
named 'workdir' file_save directory_save
aprepro
responses
objective_functions = 2
no_gradients
no_hessians
計算時間は上より圧倒的に早い。
出来たfinaldata?.datがパレート解と思われるが、その他にもいくつか出力されるので、これらのdatファイルの意味を、今後確認する必要がある。また、パレート解の評価方法についても勉強するべきだろう。
