Simple Optimization using Simultaneous DerivativesΒΆ

Consider a set of points in 2-d space that are to be arranged along a circle such that the radius of the circle is maximized, subject to constraints. We start out with the points randomly distributed within a unit circle centered about the origin. The locations of our points are determined by the values of the x and y arrays defined in our problem.

from openmdao.api import Problem, IndepVarComp, ExecComp, ScipyOptimizeDriver
import numpy as np

# note: size must be an even number
SIZE = 10
p = Problem()

indeps = p.model.add_subsystem('indeps', IndepVarComp(), promotes_outputs=['*'])

# the following were randomly generated using np.random.random(10)*2-1 to randomly
# disperse them within a unit circle centered at the origin.
indeps.add_output('x', np.array([ 0.55994437, -0.95923447,  0.21798656, -0.02158783,  0.62183717,
                                  0.04007379,  0.46044942, -0.10129622,  0.27720413, -0.37107886]))
indeps.add_output('y', np.array([ 0.52577864,  0.30894559,  0.8420792 ,  0.35039912, -0.67290778,
                                  -0.86236787, -0.97500023,  0.47739414,  0.51174103,  0.10052582]))
indeps.add_output('r', .7)

p.model.add_subsystem('circle', ExecComp('area=pi*r**2'))

p.model.add_subsystem('r_con', ExecComp('g=x**2 + y**2 - r',
                                        g=np.ones(SIZE), x=np.ones(SIZE), y=np.ones(SIZE)))

thetas = np.linspace(0, np.pi/4, SIZE)
p.model.add_subsystem('theta_con', ExecComp('g=arctan(y/x) - theta',
                                            g=np.ones(SIZE), x=np.ones(SIZE), y=np.ones(SIZE),
                                            theta=thetas))
p.model.add_subsystem('delta_theta_con', ExecComp('g = arctan(y/x)[::2]-arctan(y/x)[1::2]',
                                                  g=np.ones(SIZE//2), x=np.ones(SIZE),
                                                  y=np.ones(SIZE)))

thetas = np.linspace(0, np.pi/4, SIZE)

p.model.add_subsystem('l_conx', ExecComp('g=x-1', g=np.ones(SIZE), x=np.ones(SIZE)))

p.model.connect('r', ('circle.r', 'r_con.r'))
p.model.connect('x', ['r_con.x', 'theta_con.x', 'delta_theta_con.x'])

p.model.connect('x', 'l_conx.x')

p.model.connect('y', ['r_con.y', 'theta_con.y', 'delta_theta_con.y'])

p.driver = ScipyOptimizeDriver()
p.driver.options['optimizer'] = 'SLSQP'
p.driver.options['disp'] = False

p.model.add_design_var('x')
p.model.add_design_var('y')
p.model.add_design_var('r', lower=.5, upper=10)

# nonlinear constraints
p.model.add_constraint('r_con.g', equals=0)

IND = np.arange(SIZE, dtype=int)
ODD_IND = IND[0::2]  # all odd indices
p.model.add_constraint('theta_con.g', lower=-1e-5, upper=1e-5, indices=ODD_IND)
p.model.add_constraint('delta_theta_con.g', lower=-1e-5, upper=1e-5)

# this constrains x[0] to be 1 (see definition of l_conx)
p.model.add_constraint('l_conx.g', equals=0, linear=False, indices=[0,])

# linear constraint
p.model.add_constraint('y', equals=0, indices=[0,], linear=True)

p.model.add_objective('circle.area', ref=-1)

# setup coloring
p.driver.set_simul_deriv_color(
    (
        {
            'y': [0, 1, 0, 1, 0, 1, 0, 1, 0, 1],
            'x': [0, 1, 0, 1, 0, 1, 0, 1, 0, 1]
        },
        {
            'delta_theta_con.g': {
                'y': {
                    0: ([0, 1, 2, 3, 4], [0, 2, 4, 6, 8]),
                    1: ([0, 1, 2, 3, 4], [1, 3, 5, 7, 9])},
                'x': {
                    0: ([0, 1, 2, 3, 4], [0, 2, 4, 6, 8]),
                    1: ([0, 1, 2, 3, 4], [1, 3, 5, 7, 9])}},
            'r_con.g': {
                'y': {
                    0: ([0, 2, 4, 6, 8], [0, 2, 4, 6, 8]),
                    1: ([1, 3, 5, 7, 9], [1, 3, 5, 7, 9])},
                'x': {
                    0: ([0, 2, 4, 6, 8], [0, 2, 4, 6, 8]),
                    1: ([1, 3, 5, 7, 9], [1, 3, 5, 7, 9])}},
            'l_conx.g': {
                'x': {
                    0: ([0], [0])}},
            'theta_con.g': {
                'y': {
                    0: ([0, 1, 2, 3, 4], [0, 2, 4, 6, 8])},
                'x': {
                    0: ([0, 1, 2, 3, 4], [0, 2, 4, 6, 8])
                }
            }
        }
    )
)

p.setup(mode='fwd')
p.run_driver()

print(p['circle.area'])
[ 3.14159265]

Total derivatives with respect to x and y will be solved for simultaneously based on the color of the points shown below. Note that we have two colors and our x and y design variables are of size 10. We have a third design variable, r, that is size 1. This means that if we don’t solve for derivatives simultaneously, we must perform 21 linear solves (10 + 10 + 1) to solve for total derivatives with respect to all of our design variables. But, since each of our x and y design variables have only two colors, we can solve for all of our total derivatives using only 5 linear solves (2 + 2 + 1). This means that using simultaneous derivatives saves us 16 linear solves each time we compute our total derivatives.

Here’s our problem at the start of the optimization:

our problem before optimization

After the optimization, all of our points lie along the unit circle. The final radius is 1.0 (which corresponds to an area of PI) because we constrained our x[0] value to be 1.0.

our problem after optimization