Optimizing the Sellar Problem

In the previous tutorials we showed you how to define the Sellar model and run it directly. Now let’s see how we can optimize the Sellar problem to minimize the objective function. Here is the mathematical problem formulation for the Sellar optimization problem:

(1)\[\begin{align} \text{min}: & \ \ \ & x_1^2 + z_2 + y_1 + e^{-y_2} \\ \text{w.r.t.}: & \ \ \ & x_1, z_1, z_2 \\ \text{subject to}: & \ \ \ & \\ & \ \ \ & 3.16 - y_1 <=0 \\ & \ \ \ & y_2 - 24.0 <=0 \end{align}\]

Remember that we built our Sellar model as follows:

import numpy as np
import openmdao.api as om

from openmdao.test_suite.components.sellar import SellarDis1, SellarDis2


class SellarMDA(om.Group):
    """
    Group containing the Sellar MDA.
    """

    def setup(self):
        cycle = self.add_subsystem('cycle', om.Group(), promotes=['*'])
        cycle.add_subsystem('d1', SellarDis1(), promotes_inputs=['x', 'z', 'y2'],
                            promotes_outputs=['y1'])
        cycle.add_subsystem('d2', SellarDis2(), promotes_inputs=['z', 'y1'],
                            promotes_outputs=['y2'])

        cycle.set_input_defaults('x', 1.0)
        cycle.set_input_defaults('z', np.array([5.0, 2.0]))

        # Nonlinear Block Gauss Seidel is a gradient free solver
        cycle.nonlinear_solver = om.NonlinearBlockGS()

        self.add_subsystem('obj_cmp', om.ExecComp('obj = x**2 + z[1] + y1 + exp(-y2)',
                                                  z=np.array([0.0, 0.0]), x=0.0),
                           promotes=['x', 'z', 'y1', 'y2', 'obj'])

        self.add_subsystem('con_cmp1', om.ExecComp('con1 = 3.16 - y1'), promotes=['con1', 'y1'])
        self.add_subsystem('con_cmp2', om.ExecComp('con2 = y2 - 24.0'), promotes=['con2', 'y2'])

All the variables we need to set up the optimization are there. So now we just need the run script to execute the optimization.

SellarMDA class definition

class SellarMDA(om.Group):
    """
    Group containing the Sellar MDA.
    """

    def setup(self):
        cycle = self.add_subsystem('cycle', om.Group(), promotes=['*'])
        cycle.add_subsystem('d1', SellarDis1(), promotes_inputs=['x', 'z', 'y2'],
                            promotes_outputs=['y1'])
        cycle.add_subsystem('d2', SellarDis2(), promotes_inputs=['z', 'y1'],
                            promotes_outputs=['y2'])

        cycle.set_input_defaults('x', 1.0)
        cycle.set_input_defaults('z', np.array([5.0, 2.0]))

        # Nonlinear Block Gauss Seidel is a gradient free solver
        cycle.nonlinear_solver = om.NonlinearBlockGS()

        self.add_subsystem('obj_cmp', om.ExecComp('obj = x**2 + z[1] + y1 + exp(-y2)',
                                                  z=np.array([0.0, 0.0]), x=0.0),
                           promotes=['x', 'z', 'y1', 'y2', 'obj'])

        self.add_subsystem('con_cmp1', om.ExecComp('con1 = 3.16 - y1'), promotes=['con1', 'y1'])
        self.add_subsystem('con_cmp2', om.ExecComp('con2 = y2 - 24.0'), promotes=['con2', 'y2'])

import openmdao.api as om
from openmdao.test_suite.components.sellar_feature import SellarMDA

prob = om.Problem()
prob.model = SellarMDA()

prob.driver = om.ScipyOptimizeDriver()
prob.driver.options['optimizer'] = 'SLSQP'
# prob.driver.options['maxiter'] = 100
prob.driver.options['tol'] = 1e-8

prob.model.add_design_var('x', lower=0, upper=10)
prob.model.add_design_var('z', lower=0, upper=10)
prob.model.add_objective('obj')
prob.model.add_constraint('con1', upper=0)
prob.model.add_constraint('con2', upper=0)

# Ask OpenMDAO to finite-difference across the model to compute the gradients for the optimizer
prob.model.approx_totals()

prob.setup()
prob.set_solver_print(level=0)

prob.run_driver()

print('minimum found at')
print(prob.get_val('x')[0])
print(prob.get_val('z'))

print('minumum objective')
print(prob.get_val('obj')[0])

Optimization terminated successfully    (Exit mode 0)
            Current function value: 3.183393951729169
            Iterations: 6
            Function evaluations: 6
            Gradient evaluations: 6
Optimization Complete
-----------------------------------
minimum found at
0.0
[1.97763888e+00 8.83056605e-15]
minumum objective
3.183393951729169

Controlling the Solver Print Output

Notice the call to prob.set_solver_print(), which sets the solver output to level 0. This is the semi-quiet setting where you will only be notified if the solver failed to converge. There are lots of ways to configure the solver print output in your model to suit your needs.

Approximate the total derivatives with finite difference

In this case we’re using the SLSQP algorithm, which is a gradient-based optimization approach. Up to this point, none of our components have provided any analytic derivatives, so we’ll just finite difference across the whole model to approximate the derivatives. This is accomplished by this line of code:

prob.model.approx_totals()

Note

We’re using finite difference here for simplicity, but for larger models, finite differencing results in a high computational cost, and can have limited accuracy. It’s much better to use analytic derivatives with your models. You can learn more about that in the Advanced User Guide.