포물면 최소화 문제

아래의 최소화 문제는 OpenMDAO를 사용하여 해결합니다.

\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를 정의합니다.

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 메서드가 없더라도 계산할 수 있으며 야코비 행렬에 해당하는 사전을 반환합니다.unknowns는 $\frac {\partial f (x, y)} {\partial x} $의 실제 계산입니다.

Problem(문제) 설정

Paraboloid는 $ f (x, y) $ 함수를 나타내는 파트 (클래스)입니다.
최적화를 위해서는 설계 변수와 목적 함수를 최적화 방법을 결정해야합니다.
paraboloid.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']))

우선 1~4행째로 필요한 클래스를 임포트하고 있다
여섯 번째 줄은 top이라는 문제 (Problem의 인스턴스)를 정의합니다. 문자 그대로의이 최적화 문제의 top입니다.
7 번째 줄에서 top이라는 문제의 루트에 새로운 그룹을 만들었습니다.
9,10 행에 변수 linearize 또는 J['f_xy', 'x']가있는 p1, p2의 Component (IndepVarComp)가 추가되었습니다.
이 x는 자유 값을 취할 수있는이 문제의 디자인 변수입니다.
11 번째 줄에서 p라는 Paraboloid 인스턴스를 루트에 추가했습니다.
라인 12,13에서 설계 변수 y가 Paraboloid 입력 변수 p1.x,p2.y에 연결됩니다.
이렇게하면 설계 변수가 p1.x,p2.y를 변경하면 Paraboloid 출력 변수 p.x,p.y도 변경됩니다.

15 행 p1.x,p2.y부터는 최적화 기법을 정의합니다.
15 행에서 최적화 드라이버에 ScipyOptimizer를 지정하고 scipy로 구현 된 각 최적화 기술을 사용할 수 있습니다.
16 행째에 최적화 방법을 SLSQP. 순차 2 차 계획법이다.
라인 17,18은 설계 변수 p.f_xy의 가능한 범위를 설정합니다.
그런 다음 19 행에서 목적 함수에 araboloid 출력 변수 top.driver = ScipyOptimizer()를 지정합니다.

라인 21-24는 드라이버의 작동을 기록하는 레코더 설정입니다.
21 행 SqliteRecorder의 인수는 기록 된 데이터의 저장 파일입니다.

26 행에서 문제 해결을 27 행에서 최적화를 수행하고 28 행에서 정리했습니다.

최적화 결과

다음과 같은 실행 결과가 표시됩니다.

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)을 읽습니다.
인터프리터에서 다음 코드를 실행하십시오.

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"]

Ipython 계속

plt.plot(np.arange(0,5,1),a,"-o")
plt.xlabel("iterations")
plt.ylabel("f_xy")
plt.show()