Kriging応答曲面
理論的にはよく理解してはいないが,散らばりのあるサンプル点から曲面を推定する手法である.
OpenMDAOでラテン超方格サンプリングに示したラテン超方格サンプリングと相性が良い(散らばりのあるサンプリング)ため,工学的にはCAEを利用して応答曲面を推定する際しばしば利用される.
OpenMDAOで放物面のKriging応答曲面
取組みの概略は次のとおりである。まずOpenMDAOでラテン超方格サンプリングで実験したデータ（doe_paraboloidファイル）を読み込む.
次に読み込んだ実験データを用いてKriging近似モデルを鍛え,作成する.
最後に,作成したKriging近似モデルによる内挿値を理論解(ワイヤーフレーム)プロットすることで,モデルの妥当性を確認する.尚,次式が理論解である.
\begin{align}
& f(x,y) = (x-3)^2 + xy + (y+4)^2 - 3 \\
{\rm subject \: to} \: \: \:& -50.0\leq x \leq 50.0 \\
& -50.0\leq y \leq 50.0
\end{align}
#Groupの準備
下記のtraining_mm.pyファイルを準備する
from __future__ import print_function
from openmdao.api import Group, MetaModel, FloatKrigingSurrogate
class TrainingMM(Group):
''' FloatKriging gives responses as floats '''
def __init__(self):
super(TrainingMM, self).__init__()
# Create meta_model for f_x as the response
mm = self.add("parabo_mm", MetaModel())
mm.add_param('x', val=0.)
mm.add_param('y', val=0.)
mm.add_output('f_xy', val=0., surrogate=FloatKrigingSurrogate())
MetaModelというの近似モデルの訓練と値の評価を行うComponentをを追加している.
f_xyを近似するためにFloatKrigingSurrogate というSurrogateModelをセットしている.
FloatKrigingSurrogate単体でも近似モデルの訓練と値の評価を行う機能を有しているが,
MetaModelを用いれば複数の近似モデルを同時に訓練できると、公式ドキュメントに書かれている.
#Problem(問題）の設定
下記kriging_paraboloid.pyを準備する.前半は実験データの読み込み～Kriging近似モデルの訓練までの処理である.
#! /bin/python
import pickle
import numpy as np
import sqlitedict
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from openmdao.api import Problem
from training_mm import TrainingMM
db =sqlitedict.SqliteDict("doe_paraboloid","iterations")
res = np.array([[0.,0.,0.]] * len(db.items()))
for i, v in enumerate(db.items()):
res[i,0] = v[1]["Unknowns"]["x"]
res[i,1] = v[1]["Unknowns"]["y"]
res[i,2] = v[1]["Unknowns"]["f_xy"]
prob = Problem()
prob.root = TrainingMM()
prob.setup()
prob["parabo_mm.train:x"] = res[:,0]
prob["parabo_mm.train:y"] = res[:,1]
prob["parabo_mm.train:f_xy"] = res[:,2]
prob.run()
後半部を下記に示す。作成したkriging近似モデルと理論解の比較をプロットにて行っている
x = np.arange(-50., 50., 4.)
y = np.arange(-50., 50., 4.)
sg = prob.root.parabo_mm._init_unknowns_dict["f_xy"]["surrogate"]
f = open("./kriging_parabo","w")
pickle.dump(sg, f)
f.close()
f = open("./kriging_parabo_mm","w")
pickle.dump(prob.root.parabo_mm, f)
f.close()
xyzkrig = np.array([[xi,yi,sg.predict(np.array([xi,yi]))[0]] \
for xi in x for yi in y])
#xyzkrig = np.array([[0.,0.,0.]]*(25*25))
#cnt = 0
#for xi in x:
# for yi in y:
# xyzkrig[cnt,0] = prob["parabo_mm.x"] = xi
# xyzkrig[cnt,1] = prob["parabo_mm.y"] = yi
# prob.run()
# xyzkrig[cnt,2] = prob["parabo_mm.f_xy"]
# cnt += 1
fig = plt.figure(figsize=(6,4)); ax = Axes3D(fig)
X, Y = np.meshgrid(x, y)
Z = (X-3.0)**2. + X*Y + (Y+4.0)**2. - 3.0
ax.plot_wireframe(X,Y,Z,label="theoretical surface")
ax.scatter3D(xyzkrig[:,0], xyzkrig[:,1], xyzkrig[:,2], c="r", marker="o",label="Interpolated value")
ax.set_xlabel('x'); ax.set_ylabel('y'); ax.set_zlabel('f_xy')
plt.legend()
plt.show()
kriging_paraboloid.pyの前半部は解説は省略し後半部のみ解説する.
まず冒頭ではKriging近似モデルを評価するサンプル点を作成している.
これはLHS(ラテン超方格サンプリング)した点とは別の点であり,
$-50 \leq x,y \leq 50$の範囲を25分割した25x25のサンプル点にて評価している.
sg = prob.root.parabo_mm._init_unknowns_dict["f_xy"]["surrogate"]でプライベートクラス変数(_init_unknowns_dict)にアクセスして,MeataModelに設置されているsurrogateモデルを読み込んでいる.
褒められた処置ではないがこれには２点理由がある.
- 鍛えた近似モデルの保存が実装されていない.今回はpickeleでシリアライズした.
- MetaModelで近似モデルの値を評価しようとすると少し煩雑になる(
#xyzkrig = ・・・以下のコメント行）.
残りはプロット処理である.
理論解をワイヤーフレームにKriging近似モデルによる評価値をドットでプロットしている.
Kriging近似モデルによる値の内挿
LHSのサンプル数が多すぎたこともあり,かなり精度はよさそう.
