OpenMDAOで放物面の最小化問題を解く

# 放物面の最小化問題
下記の最小化問題をOpenMDAOを使って解く.

```math
\begin{align}
    {\rm min} \: \: \:& 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 \\
\\
    {\rm answer} \: \: \: &  f(x,y)=-27.333  \: \:  {\rm at}\:   x=6.667, \: y=-7.333 \\
\\
\end{align}
```

---
# Componentの準備
下記のようなComponent Classを継承したParaboloid Classを定義する.

```py:paraboloid.py
from openmdao.api import Component

class Paraboloid(Component):
    """ Evaluates the equation f(x,y) = (x-3)^2 + xy + (y+4)^2 - 3 """
    def __init__(self):
        super(Paraboloid, self).__init__()
        self.add_param('x', val=0.0)
        self.add_param('y', val=0.0)
        self.add_output('f_xy', shape=1)

    def solve_nonlinear(self, params, unknowns, resids):
        """f(x,y) = (x-3)^2 + xy + (y+4)^2 - 3
        """
        x = params['x']; y = params['y']
        unknowns['f_xy'] = (x-3.0)**2 + x*y + (y+4.0)**2 - 3.0

    def linearize(self, params, unknowns, resids):
        """ Jacobian for our paraboloid."""
        x = params['x']; y = params['y']
        J = {}
        J['f_xy', 'x'] = 2.0*x - 6.0 + y
        J['f_xy', 'y'] = 2.0*y + 8.0 + x
        return J
```

---
`__init__`メソッドで, 初期値0.0の`x,y`の入力変数を追加. 
shape=1 （数値型）の未知の出力変数`f_xy`を定義.

`solve_nonlinear`メソッドで$f(x,y)$の計算を行っている.ただし計算に使用する`x,y`は,引数のparamsという辞書で引き渡された値を使い,同じく引数の`unknowns`辞書の中身を更新している

`linearize`メソッドは無くても計算可能.ヤコビ行列に相当する辞書を返している.
`J['f_xy', 'x']`は$\frac{\partial  f(x,y)}{\partial  x}$の実際の計算である.

# Problem（問題）の設定
Paraboloidは$f(x,y)$関数を表す部品(クラス)である.
最適化を行うためには,設計変数や目的関数を最適化手法を決める必要がある.
paraboloid.pyを保存したディレクトリで下記のコードをスクリプトもしくはインタプリタで実行する

---

```py:opt_paraboloid.py
from __future__ import print_function
from openmdao.api import IndepVarComp, Component, Problem, Group, SqliteRecorder
from openmdao.api import ScipyOptimizer
from paraboloid import Paraboloid

top = Problem()
root = top.root = Group()

root.add('p1', IndepVarComp('x', 13.0))
root.add('p2', IndepVarComp('y', -14.0))
root.add('p', Paraboloid())
root.connect('p1.x', 'p.x')
root.connect('p2.y', 'p.y')

top.driver = ScipyOptimizer()
top.driver.options['optimizer'] = 'SLSQP'
top.driver.add_desvar('p1.x', lower=-50, upper=50)
top.driver.add_desvar('p2.y', lower=-50, upper=50)
top.driver.add_objective('p.f_xy')

recorder = SqliteRecorder('paraboloid')
recorder.options['record_params'] = True
recorder.options['record_metadata'] = True
top.driver.add_recorder(recorder)

top.setup()
top.run()
top.cleanup() 

print('Minimum of %f found at (%f, %f)' % (top['p.f_xy'], top['p.x'], top['p.y']))

```

---
まず１～４行目で必要なクラスをインポートしている
6行目でtopという問題(Problemのインスタンス)を定義している.文字通りの今回の最適化問題のtopである.
7行目でtopという問題のrootに新規グループを作っている.
9,10行目で変数`x`や`y`を持つp1,p2のComponent(IndepVarComp)を追加している.
この`p1.x,p2.y`が自由な値を取りうる本問題(top)の設計変数となる.
11行目で,pという名前のParaboloidインスタンスをrootに追加している.
12,13行目で設計変数`p1.x,p2.y`をParaboloidの入力変数`p.x,p.y`に接続している.
これにより設計変数を`p1.x,p2.y`を変化させれば,Paraboloidの出力変数`p.f_xy`も変化する.

15行目`top.driver = ScipyOptimizer()`以降では最適化手法について定義している.
15行目で最適化ドライバにScipyOptimizerを指定している.scipyで実装されている各最適化手法が使用できる.
16行目で最適化手法をSLSQP.逐次二次計画法である.
17,18行目では設計変数`p1.x,p2.y`の取りうる範囲を設定している.
そして19行目で目的関数にaraboloidの出力変数`p.f_xy`を指定してる.

21～24行目は、driverの運転を記録するrecorderの設定である.
21行目SqliteRecorderの引数は記録したデータの保存ファイルである.

26行目で問題(top Problem)のsetupを,27行目で最適化を実行し,28行目で後片付けをしている.

---
#最適化結果
下記のような実行結果が表示される.

```bash:stdout
Optimization terminated successfully.    (Exit mode 0)
            Current function value: [-27.33333329]
            Iterations: 4
            Function evaluations: 5
            Gradient evaluations: 4
Optimization Complete
Minimum of -27.333333 found at (6.666856, -7.333543)
```

次にSqliteRecorderで記録した運転記録(ファイル名:paraboloid)を読み込む.
インタプリタ上で以下のコードを実行.

```py:IPython
import numpy as np
from matplotlib import pyplot as plt
import sqlitedict

db =sqlitedict.SqliteDict("paraboloid","iterations")
a = np.zeros(5)
for i in range(0,5):
    a[i] = db[db.keys()[i]]["Unknowns"]["p.f_xy"]
```

---

```py:Ipythonつづき
plt.plot(np.arange(0,5,1),a,"-o")
plt.xlabel("iterations")
plt.ylabel("f_xy")
plt.show()
```
![figure_1.png](https://qiita-image-store.s3.amazonaws.com/0/132571/cbf71a63-f679-872a-3e3d-7fa67dcc3554.png)