LoginSignup
6

More than 5 years have passed since last update.

OpenMDAOでKriging応答曲面

Last updated at Posted at 2016-11-30
1 / 7

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ファイルを準備する

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近似モデルの訓練までの処理である.

kriging_paraboloid.py
#! /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近似モデルと理論解の比較をプロットにて行っている

kriging_paraboloid.py続き
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モデルを読み込んでいる.
褒められた処置ではないがこれには2点理由がある.
* 鍛えた近似モデルの保存が実装されていない.今回はpickeleでシリアライズした.
* MetaModelで近似モデルの値を評価しようとすると少し煩雑になる(#xyzkrig = ・・・以下のコメント行).

残りはプロット処理である.
理論解をワイヤーフレームにKriging近似モデルによる評価値をドットでプロットしている.

Kriging近似モデルによる値の内挿

LHSのサンプル数が多すぎたこともあり,かなり精度はよさそう.


kriging_paraboloid.png

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
6