Tutorial: MDAO Architectures

This tutorial shows how to create a model to solve a simple problem consisting of two coupled disciplines using several MDAO strategies, including:

  1. Multidisciplinary Design Feasible (MDF)
  2. Independent Design Feasible (IDF)
  3. Collaborative Optimization (CO)

The tutorial will introduce you to some new topics that include using the iteration hierarchy to set up models with nested optimization loops, using a solver to “close the loop” in a coupled multidisciplinary simulation, and using a broadcaster to set the values of design variables in multiple places at one time.

This tutorial covers some of the more advanced capabilities of OpenMDAO. You should read and understand the simple tutorial problem before starting this one. An understanding of the material presented in A More Complex Tutorial Problem is also recommended.

All of these tutorials use the Sellar problem which consists of two disciplines as follows:

../_images/SellarResized.png

Variables z1, z2, and x1 are the design variables over which we’d like to minimize the objective. Both disciplines are functions of z1 and z2, so they are called the global design variables, while only the first discipline is a function of x1, so it is called the local design variable. The two disciplines are coupled by the coupling variables y1 and y2. Discipline 1 takes y2 as an input, and computes y1 as an output, while Discipline 2 takes y1 as an input and computes y2 as an output. As such, the two disciplines depend on each other’s output, so iteration is required to find a set of coupling variables that satisfies both equations.

Disciplines 1 and 2 were implemented in OpenMDAO as components.

from openmdao.main.api import Component
from openmdao.lib.datatypes.api import Float

class SellarDiscipline1(Component):
    """Component containing Discipline 1"""

    # pylint: disable-msg=E1101
    z1 = Float(0.0, iotype='in', desc='Global Design Variable')
    z2 = Float(0.0, iotype='in', desc='Global Design Variable')
    x1 = Float(0.0, iotype='in', desc='Local Design Variable')
    y2 = Float(0.0, iotype='in', desc='Disciplinary Coupling')

    y1 = Float(iotype='out', desc='Output of this Discipline')


    def execute(self):
        """Evaluates the equation
        y1 = z1**2 + z2 + x1 - 0.2*y2"""

        z1 = self.z1
        z2 = self.z2
        x1 = self.x1
        y2 = self.y2

        self.y1 = z1**2 + z2 + x1 - 0.2*y2


class SellarDiscipline2(Component):
    """Component containing Discipline 2"""

    # pylint: disable-msg=E1101
    z1 = Float(0.0, iotype='in', desc='Global Design Variable')
    z2 = Float(0.0, iotype='in', desc='Global Design Variable')
    y1 = Float(0.0, iotype='in', desc='Disciplinary Coupling')

    y2 = Float(iotype='out', desc='Output of this Discipline')


    def execute(self):
        """Evaluates the equation
        y1 = y1**(.5) + z1 + z2"""

        z1 = self.z1
        z2 = self.z2

        # Note: this may cause some issues. However, y1 is constrained to be
        # above 3.16, so lets just let it converge, and the optimizer will
        # throw it out
        y1 = abs(self.y1)

        self.y2 = y1**(.5) + z1 + z2

SellarDiscipline2 contains a square root of variable y1 in its calculation. For negative values of y1, the result would be imaginary, so the absolute value is taken before the square root is applied. This component is clearly not valid for y1 < 0, and our first thought was to add a low attribute to the variable definition for y1. However, the solver that was used to converge the two disciplines occasionally forced y1 to go slightly negative. The inclusion of the absolute value solved the problem without impacting the eventual convergence of the solver.

These two components are contained in the file disciplines.py.

Reference:

Sellar, R. S., Batill, S. M., and Renaud, J. E., “Response Surface Based, Concurrent Subspace Optimization for Multidisciplinary System Design,” Proceedings References 79 of the 34th AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, January 1996.

Multidisciplinary Design Feasible (MDF)

In a Multidisciplinary Design Feasible (MDF) problem, the disciplines are directly coupled via some kind of solver, and the design variables are optimized in a single loop. The following diagram illustrates the data flow for MDF applied to the Sellar problem.

../_images/Arch-MDF.png

Data Flow for MDF Applied to the Sellar Prroblem

This diagram introduces a component called a Broadcaster. A Broadcaster is a component that enables a design variable to be set to the same value at multiple locations. If you recall, a driver such as the CONMIN optimizer contains a list of parameters, where each parameter is a location in OpenMDAO’s data hierarchy. Each parameter is a single design variable, and there is no way to indicate that one design variable might be needed at multiple component inputs in the model. We can overcome this by creating a component that passes an input to its output. Thus, CONMIN can set the design variable in this Broadcaster, and when the Broadcaster executes, the new value gets passed to all of the components that need it.

OpenMDAO doesn’t have a built-in Broadcaster, so we need to make our own. It’s a simple component with some inputs, some outputs, and an execute function that passes the inputs to the outputs.

from openmdao.main.api import Component
from openmdao.lib.datatypes.api import Float

class Broadcaster(Component):
    """Component that holds some design variables.
    This is only needed because we can't hook an optimizer up to multiple
    locations of the same design variable"""

    # pylint: disable-msg=E1101
    z1_in = Float(0.0, iotype='in', desc='Global Design Variable')
    z2_in = Float(0.0, iotype='in', desc='Global Design Variable')
    x1_in = Float(0.0, iotype='in', desc='Local Design Variable for CO')
    y1_in = Float(0.0, iotype='in', desc='Coupling Variable')
    y2_in = Float(0.0, iotype='in', desc='Coupling Variable')
    z1 = Float(0.0, iotype='out', desc='Global Design Variable')
    z2 = Float(0.0, iotype='out', desc='Global Design Variable')
    x1 = Float(0.0, iotype='out', desc='Local Design Variable for CO')
    y1 = Float(0.0, iotype='out', desc='Coupling Variable')
    y2 = Float(0.0, iotype='out', desc='Coupling Variable')

    def execute(self):
        """ Pass everything through"""
        self.z1 = self.z1_in
        self.z2 = self.z2_in
        self.x1 = self.x1_in
        self.y1 = self.y1_in
        self.y2 = self.y2_in

We’ve added the coupling variables in our Broadcaster as well, foreseeing the need for them in some of the other MDAO architectures.

The diagram also shows a solver that takes the output of the component dataflow and feeds it back into the input. OpenMDAO presently has two solvers: FixedPointIterator and BroydenSolver. The FixedPointIterator is a solver that performs fixed point iteration, which means that it keeps driving x_new = f(x_old) until convergence is achieved. In other words, y2 is passed from the output of SellarDiscipline2 to the input of SellarDiscipline1, and the loop keeps executing until the change in the value of y2 between iterations is smaller than a tolerance. The BroydenSolver is a solver based on a quasi-Newton-Raphson algorithm that uses a Broyden update to approximate the Jacobian. This solver reads the output and calculates a new input each iteration. Convergence is achieved when the residual between the output and input is driven to zero.

The major difference between the MDF problem and previous examples is the presence of nested drivers. Drivers can be nested in OpenMDAO using WorkFlows in the iteration hierarchy. A WorkFlow is an object that determines execution order for a group of Components. Each driver contains a single WorkFlow. For each iteration, a Driver will execute one pass through the WorkFlow, executing the components contained therein in the order the WorkFlow prescribes. Although in many cases a WorkFlow contains just Components, it can also contain Drivers. This allows nested iterative processes to be created. The following diagram shows an iteration hierarchy for the MDF problem.

../_images/Arch-MDF-OpenMDAO.png

Iteration Hierarchy for the MDF Problem

In the top left of this diagram, the gray box labeled Optimizer is the top level (or outermost) driver. This driver has a workflow that contains two objects – the Broadcaster and a Solver – so each time the optimizer runs an iteration, both of these components run. The Solver also has a workflow which contains the two discipline components. With the nesting of the drivers we get the behavior we want, namely, that for each optimizer iteration, the solver runs the discipline components until they converge. We now have a nested driver loop.

The execution order is determined by the components’ dataflow. Here, the broadcaster feeds the design variables to the discipline components, which are contained in the solver’s workflow, so the broadcaster must run first. Also, the data connection between the two discipline components means that SellarDiscipline1 runs before SellarDiscipline2. Sometimes a workflow may contain components that are not directly connected and can be run concurrently. Future tutorials will demonstrate this.

Now, let’s create the assembly for the MDF problem. First, we’ll create the top level optimization loop.

from openmdao.examples.mdao.disciplines import SellarDiscipline1, \
                                               SellarDiscipline2
from openmdao.examples.mdao.broadcaster import Broadcaster

from openmdao.main.api import Assembly, set_as_top
from openmdao.lib.drivers.api import CONMINdriver, FixedPointIterator

class SellarMDF(Assembly):
    """ Optimization of the Sellar problem using MDF
    Disciplines coupled with FixedPointIterator.
    """

    def __init__(self):
        """ Creates a new Assembly with this problem

        Optimal Design at (1.9776, 0, 0)

        Optimal Objective = 3.18339"""

        # pylint: disable-msg=E1101
        super(SellarMDF, self).__init__()

        # create Optimizer instance
        self.add('driver', CONMINdriver())

        # Outer Loop - Global Optimization
        self.add('bcastr', Broadcaster())
        self.add('fixed_point_iterator', FixedPointIterator())
        self.driver.workflow.add(['bcastr', 'fixed_point_iterator'])

So far nothing is really new in terms of syntax. Note that the top level driver is always named ‘driver’. However, all other drivers can be given any valid name. For this model, we’ve chosen to use the FixedPointIterator.

Next, we need to create the workflow for the solver. We create instances of SellarDiscipline1 and SellarDiscipline2 and add them to the assembly. Then, instead of adding them to the workflow of 'driver', we add them to the workflow of 'fixed_point_iterator'.

# Inner Loop - Full Multidisciplinary Solve via fixed point iteration
self.add('dis1', SellarDiscipline1())
self.add('dis2', SellarDiscipline2())
self.fixed_point_iterator.workflow.add(['dis1', 'dis2'])

Now the iteration hierarchy is finished. We still need to hook up the data connections and set up the CONMIN optimization and the fixed point iteration.

We need one connection between 'dis1' and 'dis2'. We also need to hook up 'bcastr' so that the design variables carry through to the discipline components.

self.connect('bcastr.z1','dis1.z1')
self.connect('bcastr.z1','dis2.z1')
self.connect('bcastr.z2','dis1.z2')
self.connect('bcastr.z2','dis2.z2')
self.connect('dis1.y1','dis2.y1')

Next, the parameters for the fixed point iterator must be set. FixedPointIterator is a specialized solver that is applicable only to single-input/single-output problems. As such, it does not conform to the standard driver interface. The output from SellarDiscipline2 is 'dis2.y2'. During iteration, this is the variable that is going to be sent to the input of SellarDiscipline1, which is 'dis1y2'. The parameter x_out takes the output variable while the parameter x_in takes the input variable. These are expression strings, but fixed point iteration doesn’t make sense using anything other than a single input and output. We also set the maximum number of iterations and a convergence tolerance.

# Iteration loop
self.fixed_point_iterator.x_out = 'dis2.y2'
self.fixed_point_iterator.x_in = 'dis1.y2'
self.fixed_point_iterator.max_iteration = 1000
self.fixed_point_iterator.tolerance = .0001

Finally, the CONIM optimization is set up.

# Optimization parameters
self.driver.add_objective('(dis1.x1)**2 + bcastr.z2 + dis1.y1 + math.exp(-dis2.y2)')

self.driver.add_parameter('bcastr.z1_in', low = -10.0, high = 10.0)
self.driver.add_parameter('bcastr.z2_in', low = 0.0,   high = 10.0)
self.driver.add_parameter('dis1.x1',      low = 0.0,   high = 10.0)

self.driver.add_constraint('3.16 < dis1.y1')
self.driver.add_constraint('dis2.y2 < 24.0')

self.driver.cons_is_linear = [1, 1]
self.driver.iprint = 0
self.driver.itmax = 30
self.driver.fdch = .001
self.driver.fdchm = .001
self.driver.delfun = .0001
self.driver.dabfun = .000001
self.driver.ctlmin = 0.0001

The process of getting the optimizer to reach a value close enough to the correct optimum required the use of more of CONMIN’s settings. The fdchm parameter is the minimum absolute step size that the finite difference uses, and fdch is the step size relative to the design variable. Dabfun is the absolute change in the objective function to indicate convergence (i.e., if the objective function changes by less than dabfun, then the problem is converged). Similarly, delfun is the relative change of the objective function with respect to the value at the previous step. Finally, ctlmin is the minimum constraint thickness for the linear constraints. We also use cons_is_linear to let CONMIN know that both constraints are linear. This can speed up the algorithm, though it hardly matters here.

As before, the add_constraint method is used to add our constraints. This time however, we used a more general expression for the first constraint. Expression strings in OpenMDAO can also be parsed as inequalities, so all of the following are equivalent ways of defining this constraint:

self.driver.add_constraint('3.16 - dis1.y1 < 0')
self.driver.add_constraint('3.16 < dis1.y1')
self.driver.add_constraint('-3.16 > -dis1.y1')

Finally, putting it all together gives:

from openmdao.examples.mdao.disciplines import SellarDiscipline1, \
                                               SellarDiscipline2
from openmdao.examples.mdao.broadcaster import Broadcaster

from openmdao.main.api import Assembly, set_as_top
from openmdao.lib.drivers.api import CONMINdriver, FixedPointIterator

class SellarMDF(Assembly):
    """ Optimization of the Sellar problem using MDF
    Disciplines coupled with FixedPointIterator.
    """

    def __init__(self):
        """ Creates a new Assembly with this problem

        Optimal Design at (1.9776, 0, 0)

        Optimal Objective = 3.18339"""

        # pylint: disable-msg=E1101
        super(SellarMDF, self).__init__()

        # create Optimizer instance
        self.add('driver', CONMINdriver())

        # Outer Loop - Global Optimization
        self.add('bcastr', Broadcaster())
        self.add('fixed_point_iterator', FixedPointIterator())
        self.driver.workflow.add(['bcastr', 'fixed_point_iterator'])

        # Inner Loop - Full Multidisciplinary Solve via fixed point iteration
        self.add('dis1', SellarDiscipline1())
        self.add('dis2', SellarDiscipline2())
        self.fixed_point_iterator.workflow.add(['dis1', 'dis2'])

        # Make all connections
        self.connect('bcastr.z1','dis1.z1')
        self.connect('bcastr.z1','dis2.z1')
        self.connect('bcastr.z2','dis1.z2')
        self.connect('bcastr.z2','dis2.z2')
        self.connect('dis1.y1','dis2.y1')

        # Iteration loop
        self.fixed_point_iterator.x_out = 'dis2.y2'
        self.fixed_point_iterator.x_in = 'dis1.y2'
        self.fixed_point_iterator.max_iteration = 1000
        self.fixed_point_iterator.tolerance = .0001

        # Optimization parameters
        self.driver.add_objective('(dis1.x1)**2 + bcastr.z2 + dis1.y1 + math.exp(-dis2.y2)')

        self.driver.add_parameter('bcastr.z1_in', low = -10.0, high = 10.0)
        self.driver.add_parameter('bcastr.z2_in', low = 0.0,   high = 10.0)
        self.driver.add_parameter('dis1.x1',      low = 0.0,   high = 10.0)

        self.driver.add_constraint('3.16 < dis1.y1')
        self.driver.add_constraint('dis2.y2 < 24.0')

        self.driver.cons_is_linear = [1, 1]
        self.driver.iprint = 0
        self.driver.itmax = 30
        self.driver.fdch = .001
        self.driver.fdchm = .001
        self.driver.delfun = .0001
        self.driver.dabfun = .000001
        self.driver.ctlmin = 0.0001

This problem is contained in sellar_MDF.py. Executing it at the command line should produce output that resembles this:

$ python sellar_MDF.py
CONMIN Iterations:  12
Minimum found at (1.977657, 0.000000, 0.000000)
Couping vars: 3.160068, 3.755315
Minimum objective:  3.18346116811
Elapsed time:  0.121051073074 seconds

We can also replace the fixed point iterator with a better solver. Fixed point iteration works for some problems, including this one, but may not converge to a solution for other problems. OpenMDAO also contains a Broyden solver called BroydenSolver. This solver is based on a quasi-Newton-Raphson algorithm found in scipy.nonlinear. It uses a Broyden update to approximate the Jacobian. If we replace FixedPointIterator with BroydenSolver, the optimizer’s workflow looks like this:

# Don't forget to put the import in your header
from openmdao.lib.drivers.api import BroydenSolver

# Outer Loop - Global Optimization
self.add('bcastr', Broadcaster())
self.add('solver', BroydenSolver())
self.driver.workflow.add(['bcastr', 'solver'])

Next, we set up our parameters for the inner loop. The Broyden solver can be connected using the standard driver interface. It can take multiple inputs and outputs though we only have one input and one output in this example.

# Iteration loop
self.solver.add_parameter('dis1.y2', low=-9.e99, high=9.e99)
self.solver.add_constraint('dis2.y2 = dis1.y2')
self.solver.itmax = 10
self.solver.alpha = .4
self.solver.tol = .0000001
self.solver.algorithm = "broyden2"

The input is selected using add_parameter. You might also be familiar with the term independent used to describe this. Here, we’ve given a low and a high attribute, but we’ve set them to very large negative and positive values as the Broyden solver doesn’t use either of these. The output is specified by adding an equality constraint. A solver essentially tries to drive something to zero. In this case, we want to drive the residual error in the coupled variable y2 to zero. An equality constraint is defined with an expression string which is parsed for the equals sign, so the following constraints are equivalent:

# Iteration loop
self.solver.add_constraint('dis2.y2 = dis1.y2')
self.solver.add_constraint('dis2.y2 - dis1.y2 = 0')

Equality constraints may also be available for some optimizers, but you should verify that they are supported. CONMIN does not support equality constraints.

Individual Design Feasible (IDF)

Next, we will look at how to set up the Individual Design Feasible (IDF) architecture using the Sellar problem. In IDF, the coupling between the disciplines is removed, and the input coupling variables are added to the optimizer’s design variables. The algorithm calls for two new equality constraints that constrain to zero the residual error between the coupling variable output by the optimizer and the coupling variable output by the components. This assures that the solution is a feasible coupling, though it is achieved through the optimizer’s additional effort instead of a solver. The data flow for IDF is illustrated in the following diagram:

../_images/Arch-IDF.png

Data Flow for IDF

IDF needs only one driver, so there is just one workflow. The broadcaster and the two disciplines are executed sequentially.

../_images/Arch-IDF-OpenMDAO.png

Iteration Hierarchy for IDF

Next, we will create the SellarIDF assembly. First, all of our components are instantiated and the workflow is defined.

from openmdao.examples.mdao.disciplines import SellarDiscipline1, \
                                               SellarDiscipline2
from openmdao.examples.mdao.broadcaster import Broadcaster

from openmdao.main.api import Assembly, set_as_top
from openmdao.lib.drivers.api import CONMINdriver


class SellarIDF(Assembly):
    """ Optimization of the Sellar problem using IDF"""

    def __init__(self):
        """ Creates a new Assembly with this problem

        Optimal Design at (1.9776, 0, 0)

        Optimal Objective = 3.18339"""

        # pylint: disable-msg=E1101

        super(SellarIDF, self).__init__()

        # create Optimizer instance
        self.add('driver', CONMINdriver())

        # Disciplines
        self.add('bcastr', Broadcaster())
        self.add('dis1', SellarDiscipline1())
        self.add('dis2', SellarDiscipline2())

        # Driver process definition
        self.driver.workflow.add(['bcastr', 'dis1', 'dis2'])

        # Make all connections
        self.connect('bcastr.z1','dis1.z1')
        self.connect('bcastr.z1','dis2.z1')
        self.connect('bcastr.z2','dis1.z2')
        self.connect('bcastr.z2','dis2.z2')

We’ve also hooked up our data connections. Only the design variables that are shared by both components need to be connected to the broadcaster.

All that is left to do is set up the CONMIN optimizer.

# Optimization parameters
self.driver.add_objective('(dis1.x1)**2 + bcastr.z2 + dis1.y1 + math.exp(-dis2.y2)')

self.driver.add_parameter('bcastr.z1_in', low = -10.0, high=10.0)
self.driver.add_parameter('bcastr.z2_in', low = 0.0,   high=10.0)
self.driver.add_parameter('dis1.x1',      low = 0.0,   high=10.0)
self.driver.add_parameter('dis2.y1',      low = 3.16,  high=10.0)
self.driver.add_parameter('dis1.y2',      low = -10.0, high=24.0)

self.driver.add_constraint('(dis2.y1-dis1.y1)**3 < 0')
self.driver.add_constraint('(dis1.y1-dis2.y1)**3 < 0')
self.driver.add_constraint('(dis2.y2-dis1.y2)**3 < 0')
self.driver.add_constraint('(dis1.y2-dis2.y2)**3 < 0')
self.driver.iprint = 0
self.driver.itmax = 100
self.driver.fdch = .003
self.driver.fdchm = .003
self.driver.delfun = .0001
self.driver.dabfun = .00001
self.driver.ct = -.01
self.driver.ctlmin = 0.001

Notice that the coupling variables are included as optimizer parameters. We also introduce the CONMIN parameter ct, which is the constraint thickness for nonlinear constraints. Our constraints are nonlinear, but note that any constraint that involves a component output is most likely a nonlinear constraint because outputs are usually nonlinear functions of the design variables.

Since CONMIN doesn’t support equality constraints, we have to fall back on a trick where we replace it with an equivalent pair of inequality constraints. For example, if we want to constrain x=2, we could constraint x<=2 and x>=2 and let the optimizer converge to a solution where both constraints are active. Stability may be questionable for such a method, so it is always advisable to use an optimizer that has equality constraints rather than trying to squeeze a solution out of an optimizer this way. In particular, be careful about trying a fancier solution such as constraining abs(dis2.y1-dis1.y1)<=0. This nonlinear constraint has a discontinuous slope, and CONMIN won’t handle that constraint very well. Here, we take (dis2.y1-dis1.y1) and turn it into a cubic expression, which seemed to make the problem a little less sensitive to changes in the computational environment (32 bit vs 64 bit, etc.)

This problem is contained in sellar_IDF.py. Executing it at the command line should produce output that resembles this:

$ python sellar_IDF.py
CONMIN Iterations:  10
Minimum found at (1.976427, 0.000287, 0.000000)
Couping vars: 3.156521, 3.754359
Minimum objective:  3.18022323743
Elapsed time:  0.200541973114 seconds

Collaborative Optimization (CO)

Next, we will set up a model that solves the Sellar problem by Collaborative Optimization (CO). CO is a two-level architecture with three optimizer loops, one at each discipline, and one acting globally. The global optimizer drives the design and coupling variables towards an optimal solution that minimizes the objective while constraining to zero the sum of the squares of the residuals between the values commanded by the global optimizer and those set by the local optimizers. Each local optimizer operates on its own discipline, driving its design variables while minimizing the residual between the actual value of the design variables and the values commanded by the global optimizer.

CO for the Sellar case is very interesting because there are no component data connections. All values are passed through the expression strings for the objectives, constraints, and parameters of the various optimizers, as shown in the next diagram.

../_images/Arch-CO.png

Data Flow for CO

The CO model has three optimizers, so there are three workflows. The top level workflow includes the broadcaster and the two lower level optimizers, and each of those optimizers has a workflow with just the discipline component. This can be seen in the next figure.

../_images/Arch-CO-OpenMDAO.png

Iteration Hierarchy for CO

First, we create the component instances and set up this iteration hierarchy. Notice that there are three drivers, and we add each component to one of the three workflows.

from openmdao.examples.mdao.disciplines import SellarDiscipline1, \
                                               SellarDiscipline2
from openmdao.examples.mdao.broadcaster import Broadcaster

from openmdao.main.api import Assembly, set_as_top
from openmdao.lib.drivers.api import CONMINdriver

class SellarCO(Assembly):
    """Solution of the sellar analytical problem using CO.
    """

    def __init__(self):
        """ Creates a new Assembly with this problem

        Optimal Design at (1.9776, 0, 0)

        Optimal Objective = 3.18339"""

        # pylint: disable-msg=E1101
        super(SellarCO, self).__init__()

        # Global Optimization
        self.add('driver', CONMINdriver())
        self.add('bcastr', Broadcaster())
        self.add('localopt1', CONMINdriver())
        self.add('localopt2', CONMINdriver())
        self.driver.workflow.add(['bcastr', 'localopt1',
                                  'localopt2'])

        # Local Optimization 1
        self.add('dis1', SellarDiscipline1())
        self.localopt1.workflow.add('dis1')

        # Local Optimization 2
        self.add('dis2', SellarDiscipline2())
        self.localopt2.workflow.add('dis2')

Notice that there are no data connections, so we never need to call self.connect.

Now we need to set up the parameters for the outer optimization loop.

#Parameters - Global Optimization
self.driver.add_objective('(bcastr.x1)**2 + bcastr.z2 + bcastr.y1' +
                                        '+ math.exp(-bcastr.y2)')
self.driver.add_parameter('bcastr.z1_in', low = -10.0, high = 10.0)
self.driver.add_parameter('bcastr.z2_in', low = 0.0,   high = 10.0)
self.driver.add_parameter('bcastr.x1_in', low = 0.0,   high = 10.0)
self.driver.add_parameter('bcastr.y1_in', low = 3.16,  high = 10.0)
self.driver.add_parameter('bcastr.y2_in', low = -10.0, high = 24.0)

con1 = '(bcastr.z1-dis1.z1)**2 + (bcastr.z2-dis1.z2)**2 + ' + \
       '(bcastr.x1-dis1.x1)**2 + ' + \
       '(bcastr.y1-dis1.y1)**2 + (bcastr.y2-dis1.y2)**2  < 0'
con2 = '(bcastr.z1-dis2.z1)**2 + (bcastr.z2-dis2.z2)**2 + ' + \
       '(bcastr.y1-dis2.y1)**2 + (bcastr.y2-dis2.y2)**2  < 0'
self.driver.add_constraint(con1)
self.driver.add_constraint(con2)

self.driver.printvars = ['dis1.y1','dis2.y2']
self.driver.iprint = 0
self.driver.itmax = 100
self.driver.fdch = .003
self.driver.fdchm = .003
self.driver.delfun = .0001
self.driver.dabfun = .00001
self.driver.ct = -.001
self.driver.ctlmin = 0.001

Here we are able to build up a complicated expression for the sum of the squares of all of the residuals and use it as our constraint. Note that this is another example of a constraint that would be better served as an equality constraint, but it’s not problematic because a sum of squares is one-sided. We have two constraints here, one for each discipline.

Finally, we set up our local optimization loops.

#Parameters - Local Optimization 1
self.localopt1.add_objective('(bcastr.z1-dis1.z1)**2 + ' + \
                           '(bcastr.z2-dis1.z2)**2 + ' + \
                           '(bcastr.x1-dis1.x1)**2 + ' + \
                           '(bcastr.y1-dis1.y1)**2 + ' + \
                           '(bcastr.y2-dis1.y2)**2')
self.localopt1.add_parameter('dis1.z1', low = -10.0, high = 10.0)
self.localopt1.add_parameter('dis1.z2', low = 0.0,   high = 10.0)
self.localopt1.add_parameter('dis1.x1', low = 0.0,   high = 10.0)
self.localopt1.add_parameter('dis1.y2', low = -10.0, high = 24.0)
self.localopt1.iprint = 0
self.localopt1.itmax = 100
self.localopt1.fdch = .003
self.localopt1.fdchm = .003
self.localopt1.delfun = .001
self.localopt1.dabfun = .00001

#Parameters - Local Optimization 2
self.localopt2.add_objective('(bcastr.z1-dis2.z1)**2 + ' + \
                             '(bcastr.z2-dis2.z2)**2 + ' + \
                             '(bcastr.y1-dis2.y1)**2 + ' + \
                             '(bcastr.y2-dis2.y2)**2')
self.localopt2.add_parameter('dis2.z1', low = -10.0, high = 10.0)
self.localopt2.add_parameter('dis2.z2', low = 0.0,   high = 10.0)
self.localopt2.add_parameter('dis2.y1', low = 3.16,  high = 10.0)
self.localopt2.iprint = 0
self.localopt2.itmax = 100
self.localopt2.fdch = .003
self.localopt2.fdchm = .003
self.localopt2.delfun = .001
self.localopt2.dabfun = .00001

This problem is contained in sellar_CO.py. Executing it at the command line should produce output that resembles this:

$ python sellar_CO.py
CONMIN Iterations:  60
Minimum found at (2.007823, 0.000000, 0.003927)
Couping vars: 3.164985, 3.802037
Minimum objective:  3.18215119048
Elapsed time:  22.5064539909 seconds

After 60 iterations of the top level optimizer, CO gives an answer that’s slightly off the optimum, but the calculated objective is still very close.