[PYTHON] Solving the paraboloid minimization problem with OpenMDAO

Parabolic minimization problem

Solve the following minimization problem using 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 preparation

Define a Paraboloid Class that inherits the Component Class as shown below.

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

In the __init __ method, add the input variables of x and y with the initial value of 0.0. Define an unknown output variable f_xy with shape = 1 (numeric type).

The solve_nonlinear method is used to calculate $ f (x, y) $. However, the x, y used in the calculation uses the value passed in the dictionary called params, which is also an argument, ʻunknowns. `Updating the contents of the dictionary

Computable without the linearize method. Returns a dictionary equivalent to the Jacobian matrix. J ['f_xy','x'] is the actual calculation of $ \ frac {\ partial f (x, y)} {\ partial x} $.

Problem settings

Paraboloid is a component (class) that represents the $ f (x, y) $ function. In order to perform optimization, it is necessary to determine the optimization method for design variables and objective functions. Execute the following code with a script or interpreter in the directory where paraboloid.py is saved

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']))

First, the required classes are imported in the 1st to 4th lines. Line 6 defines the problem (an instance of Problem), which is literally the top of this optimization problem. Line 7 creates a new group for the root in question, top. Lines 9 and 10 add a Component (IndepVarComp) of p1 and p2 with variables x and y. This p1.x, p2.y is the design variable of this problem (top) that can take any value. Line 11 adds a Paraboloid instance named p to root. Lines 12 and 13 connect the design variables p1.x, p2.y to the Paraboloid input variables p.x, p.y. As a result, if the design variables are changed to p1.x and p2.y, the output variable p.f_xy of Paraboloid will also change.

The 15th line top.driver = ScipyOptimizer () and later define the optimization method. Line 15 specifies ScipyOptimizer as the optimization driver. Each optimization method implemented in scipy can be used. The optimization method on the 16th line is SLSQP. Lines 17 and 18 set the possible range of design variables p1.x, p2.y. And on the 19th line, the output variable p.f_xy of araboloid is specified for the objective function.

Lines 21-24 are the recorder settings that record the driver's driving. Line 21 The argument of SqliteRecorder is the save file of the recorded data.

The 26th line is the setup for the top problem, the 27th line is the optimization, and the 28th line is the cleanup.

Optimization result

The following execution result is displayed.

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)

Next, read the operation record (file name: paraboloid) recorded by SqliteRecorder. Execute the following code on the interpreter.

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 continued

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