勉強会＠関西のKさんの資料に乗っていた例題を参考に、OFを利用する最適化をやってみよう（細かい寸法は勝手にアレンジしています。）
設計変数２、目的関数３、制約条件は設計変数範囲のみという基本的な問題。
簡単に説明すると、規定の入り口、出口形状が、ベンドでつながれていて、ベンドの中心線は、3次ベジエ曲線で規定されている（図中のP０～P３の4点）。このベジエ曲線の定義点に目的変数がx１、x２がふくまれている。ベンドの幅は５０mmで一定としている。このため、メッシュ生成が破綻しないように設計変数に大小の制限が加わっている（blockMesh利用）。
目的関数は、１）ベンドの長さを最小化２）入り口の平均圧力の最小化（圧力損失の最小化）３）出口の流速変化の最小化（一様性の最大化）になっている。
最適化の流れは、１）線形計画法による設計変数の選定、２）近似応答曲面の推定（surrogate model）、３）MOGAなどの手法による最適化という流れで、まずは１）から。
設計計画法だが、一番簡単に設計変数をコンスタントに区切って（以下の例ではx1, x2の範囲をそれぞれ5分割するので、6×6＝36点の設計変数ベクトルとなる）、すべての点を計算してみる。網羅的にやって、目的関数と設計変数の関係にあたりをつけるのが目的です。
なお、dakotaに入力となる、results.outというファイルには、各々のケースディレクトリ中のOFのランスクリプトで出力させてしまっている。
＜dakota_dace.in＞
# dakota -i dakota_dace.in
environment
graphics
tabular_graphics_data
tabular_graphics_file = 'multidim.dat'
method
multidim_parameter_study
partitions = 5 5
model
single
variables
continuous_design = 2
lower_bounds 6.0 6.0
upper_bounds 18.0 18.0
descriptors 'x1' "x2"
interface,
system
analysis_driver = 'run_openfoam.sh'
parameters_file = 'params.in'
results_file = 'results.out'
work_directory directory_tag
copy_files = 'pipe2D/*'
named 'workdir' file_save directory_save
aprepro
responses
objective_functions = 3
no_gradient
no_hessians
計算結果は下記の通りで、目的関数f１（ベンドの長さ）については、基本的にx1,x2とも単純な比例関係を示す（今回のベジエ曲線の定義からは、当たり前の結果と言えばその通り）。
目的関数f２、f３については、わりとスムーズな変化を示し、設計変数（８、８）くらいで、いずれも最小値が存在するような結果となった。
次に、近似応答曲面を作るのを入れずに、直接MOGAを実行してどれくらい時間がかかるかを確認してみる。
MOGAの設定は、初期値は乱数による生成、最大試行回数３０００、交配率0.8、変異率0.1、淘汰率0.9、最大４０世代とした（これらの設定値はどのような意味と、結果への影響を持つかは確認する必要がある。）
＜dakota_moga.in＞
# # dakota -i dakota_moga.in
environment
graphics
tabular_graphics_data
tabular_graphics_file = 'objective.dat'
method
moga
output silent
seed = 10983
final_solutions = 3
max_function_evaluations = 3000
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
convergence_type metric_tracker
percent_change = 0.05 num_generations = 40
model
single
variables
continuous_design = 2
initial_point 10 10
lower_bounds 6 6
upper_bounds 18 18
descriptors 'x1' 'x2'
interface
# system
fork
asynchronous
evaluation_concurrency = 4
analysis_driver = 'run_openfoam.sh'
parameters_file = 'params.in'
results_file = 'results.out'
work_directory directory_tag
copy_files = 'pipe2D/*'
named 'workdir' file_save directory_save
aprepro
responses
objective_functions = 3
no_gradients
no_hessian
計算時間は、３０００回の計算（ただし、設計変数の組を4並列で同時起動し、各OFソルバーはシリアル計算）で、約３０分かかった。（Core i7-960）
<<<<< Function evaluation summary: 3000 total (3000 new, 0 duplicate)
<<<<< Best parameters (set 1) =
7.9207336763e+00 x1
7.4609955798e+00 x2
<<<<< Best objective functions (set 1) =
3.0082169221e+01
6.1509600000e+00
2.0211800000e+03
<<<<< Best data captured at function evaluation 1136
<<<<< Best parameters (set 2) =
8.0410799133e+00 x1
7.3764506082e+00 x2
<<<<< Best objective functions (set 2) =
3.0089009353e+01
6.1500000000e+00
2.0211700000e+03
<<<<< Best data captured at function evaluation 2798
<<<<< Best parameters (set 3) =
7.7423646421e+00 x1
7.6457909186e+00 x2
<<<<< Best objective functions (set 3) =
3.0083005512e+01
6.1505300000e+00
2.0211800000e+03
<<<<< Best data captured at function evaluation 1131
<<<<< Iterator moga completed.
<<<<< Environment execution completed.
DAKOTA execution time in seconds:
Total CPU = 9.73 [parent = 9.72352, child = 0.006478]
Total wall clock = 1732.22
ということで、（精度を落とさずに）如何にこの試行計算の回数を減らすかが「肝」なようなので、近似応答曲面を用いたsurrogate-based optimizationを考えてみる。
<dakota_sbopt.in>
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 = 10
replace_points
method
id_method = 'MOGA'
moga
output silent
seed = 10983
final_solutions = 3
population_size = 500
max_function_evaluations = 3000
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
niching_type distance 0.05 0.05
postprocessor_type
orthogonal_distance 0.05 0.05
replacement_type below_limit = 6
shrinkage_percentage = 0.9
convergence_type metric_tracker
percent_change = 0.05 num_generations = 40
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 = 100 seed = 531
sample_type lhs
model
id_model = 'TRUTH'
single
variables
continuous_design = 2
initial_point 10 10
lower_bounds 6 6
upper_bounds 18 18
descriptors 'x1' 'x2'
interface
# system
fork
asynchronous
evaluation_concurrency = 4
analysis_driver = 'run_openfoam.sh'
parameters_file = 'params.in'
results_file = 'results.out'
work_directory directory_tag
copy_files = 'pipe2D/*'
named 'workdir' file_save directory_save
aprepro
responses
objective_functions = 3
no_gradients
no_hessians
計算結果は、以下の通りで非常に短時間で終了しており、近似応答曲面から、計算の試行回数も130回で済んでいる。これは、応答曲面が、非常に単調なせいもあるだろう。
ただし、結果については若干のずれがあるので、これはMOGAの設定などを詰める必要がありそう。
<<<<< Function evaluation summary (APPROX_INTERFACE): 30000 total (30000 new, 0 duplicate)
<<<<< Function evaluation summary: 130 total (130 new, 0 duplicate)
<<<<< Best parameters (set 1) =
7.8864829679e+00 x1
7.7475738105e+00 x2
<<<<< Best objective functions (set 1) =
3.0127621565e+01
6.1405100000e+00
2.0211100000e+03
<<<<< Best data captured at function evaluation 128
<<<<< Best parameters (set 2) =
8.0046206014e+00 x1
7.7475738105e+00 x2
<<<<< Best objective functions (set 2) =
3.0149252738e+01
6.1364200000e+00
2.0210800000e+03
<<<<< Best data captured at function evaluation 129
<<<<< Best parameters (set 3) =
7.7474202205e+00 x1
7.8376652830e+00 x2
<<<<< Best objective functions (set 3) =
3.0118692647e+01
6.1435200000e+00
2.0211300000e+03
<<<<< Best data captured at function evaluation 130
<<<<< Iterator surrogate_based_global completed.
<<<<< Environment execution completed.
DAKOTA execution time in seconds:
Total CPU = 61.52 [parent = 61.5126, child = 0.007352]
Total wall clock = 138.444
※ グラフのf3のz軸スケールが何かおかしい。数値的には2010とかで、gnuplotなんかは正しくプロットしてくれる。matplotlibのバグかな？
