# Paraboloid Tutorial - Simple Optimization Problem¶

This tutorial will show you how to set up a simple optimization of a paraboloid. You’ll create a paraboloid Component (with analytic derivatives), then put it into a Problem and set up an optimizer Driver to minimize an objective function.

Here is the code that defines the paraboloid and then runs it. You can copy this code into a file, and run it directly.

from __future__ import print_function

from openmdao.api import IndepVarComp, Component, Problem, Group

class Paraboloid(Component):
""" Evaluates the equation f(x,y) = (x-3)^2 + xy + (y+4)^2 - 3 """

def __init__(self):
super(Paraboloid, self).__init__()

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

if __name__ == "__main__":

top = Problem()

root = top.root = Group()

root.connect('p1.x', 'p.x')
root.connect('p2.y', 'p.y')

top.setup()
top.run()

print(top['p.f_xy'])


Now we will go through each section and explain how this code works.

## Building the component¶

from __future__ import print_function

from openmdao.api import IndepVarComp, Component, Problem, Group


We need to import some OpenMDAO classes. We also import the print_function to ensure compatibility between Python 2.x and 3.x. You don’t need the import if you are running in Python 3.x.

class Paraboloid(Component):


OpenMDAO provides a base class, Component, which you should inherit from to build your own components and wrappers for analysis codes. Components can declare three kinds of variables, parameters, outputs and states. A Component operates on its parameters to compute unknowns, which can be explicit outputs or implicit states. For the Paraboloid Component, we will only be using explicit outputs.

def __init__(self):
super(Paraboloid, self).__init__()



This code defines the input parameters of the Component, x and y, and initializes them to 0.0. These will be design variables which could be used to minimize the output when doing optimization. It also defines the explicit output, f_xy, but only gives it a shape. If shape is 1, the value is initialized to 0.0, a scalar. If shape is any other value, the value of the variable is initialized to numpy.zeros(shape, dtype=float).

def solve_nonlinear(self, params, unknowns, resids):
"""f(x,y) = (x-3)^2 + xy + (y+4)^2 - 3
Optimal solution (minimum): x = 6.6667; y = -7.3333
"""

x = params['x']
y = params['y']

unknowns['f_xy'] = (x-3.0)**2 + x*y + (y+4.0)**2 - 3.0


The solve_nonlinear method is responsible for calculating outputs for a given set of parameters. The parameters are given in the params dictionary that is passed in to this method. Similarly, the outputs are assigned values using the unknowns dictionary that is passed in.

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


The linearize method is used to compute analytic partial derivatives of the unknowns with respect to params (partial derivatives in OpenMDAO context refer to derivatives for a single component by itself). The returned value, in this case J, should be a dictionary whose keys are tuples of the form (‘unknown’, ‘param’) and whose values are n-d arrays or scalars. Just like for solve_nonlinear, the values for the parameters are accessed using dictionary arguments to the function.

The definition of the Paraboloid Component class is now complete. We will now make use of this class to run a model.

## Setting up the model¶

if __name__ == "__main__":

top = Problem()
root = top.root = Group()


An instance of an OpenMDAO Problem is always the top object for running a model. Each Problem in OpenMDAO must contain a root Group. A Group is a System that contains other Components or Groups.

This code instantiates a Problem object and sets the root to be an empty Group.

root.add('p1', IndepVarComp('x', 3.0))


Now it is time to add components to the empty group. IndepVarComp is a Component that provides the source for a variable which we can later give to a Driver as a design variable to control.

We created two IndepVarComps (one for each param on the Paraboloid component), gave them names, and added them to the root Group. The add method takes a name as the first argument, and a Component instance as the second argument. The numbers 3.0 and -4.0 are values chosen for each as starting points for the optimizer.

Note

Take care setting the initial values, as in some cases, various initial points for the optimization will lead to different results.

root.add('p', Paraboloid())


Then we add the paraboloid using the same syntax as before, giving it the name ‘p’.

root.connect('p1.x', 'p.x')
root.connect('p2.y', 'p.y')


Then we connect up the outputs of the IndepVarComps to the parameters of the Paraboloid. Notice the dotted naming convention used to refer to variables. So, for example, p1 represents the first IndepVarComp that we created to set the value of x and so we connect that to parameter x of the Paraboloid. Since the Paraboloid is named p and has a parameter x, it is referred to as p.x in the call to the connect method.

Every problem has a Driver and for most situations, we would want to set a Driver for the Problem using code like this

top.driver = SomeDriver()


For this very simple tutorial, we do not need to set a Driver, we will just use the default, built-in driver, which is Driver. ( Driver also serves as the base class for all Drivers. ) Driver is the simplest driver possible, running a Problem once.

top.setup()


Before we can run our model we need to do some setup. This is done using the setup method on the Problem. This method performs all the setup of vector storage, data transfer, etc., necessary to perform calculations. Calling setup is required before running the model.

top.run()


Now we can run the model using the run method of Problem.

print(top['p.f_xy'])


Finally, we print the output of the Paraboloid Component using the dictionary-style method of accessing variables on the problem instance. Putting it all together:

top = Problem()
root = top.root = Group()

root.connect('p1.x', 'p.x')
root.connect('p2.y', 'p.y')

top.setup()
top.run()

print(top['p.f_xy'])


The output should look like this:

-15.0


The IndepVarComp component is used to define a source for an unconnected param that we want to use as an independent variable that can be declared as a design variable for a driver. In our case, we want to optimize the Paraboloid model, finding values for ‘x’ and ‘y’ that minimize the output ‘f_xy.’

Sometimes we just want to run our component once to see the result. Similarly, sometimes we have params that will be constant through our optimization, and thus don’t need to be design variables. In either of these cases, the IndepVarComp is not required, and we can build our model while leaving those parameters unconnected. All unconnected params use their default value as the initial value. You can set the values of any unconnected params the same way as any other variables by doing the following:

top = Problem()
root = top.root = Group()

top.setup()

# Set values for x and y
top['x'] = 5.0
top['y'] = 2.0

top.run()

print(top['p.f_xy'])


This can only be done after setup is called. Note that the promoted names ‘x’ and ‘y’ are used.

The new output should look like this:

47.0


## Optimization of the Paraboloid¶

Now that we have the paraboloid model set up, let’s do a simple unconstrained optimization. Let’s find the minimum point on the Paraboloid over the variables x and y. This requires the addition of just a few more lines.

First, we need to import the optimizer.

from openmdao.api import ScipyOptimizer


The main optimizer built into OpenMDAO is a wrapper around Scipy’s minimize function. OpenMDAO supports 9 of the optimizers built into minimize. The ones that will be most frequently used are SLSQP and COBYLA, since they are the only two in the minimize package that support constraints. We will use SLSQP because it supports OpenMDAO-supplied gradients.

top = Problem()
root = top.root = Group()

# Initial value of x and y set in the IndepVarComp.

root.connect('p1.x', 'p.x')
root.connect('p2.y', 'p.y')

top.driver = ScipyOptimizer()
top.driver.options['optimizer'] = 'SLSQP'

top.setup()

# You can also specify initial values post-setup
top['p1.x'] = 3.0
top['p2.y'] = -4.0

top.run()

print('\n')
print('Minimum of %f found at (%f, %f)' % (top['p.f_xy'], top['p.x'], top['p.y']))


Every driver has an options dictionary which contains important settings for the driver. These settings tell ScipyOptimizer which optimization method to use, so here we select ‘SLSQP’. For all optimizers, you can specify a convergence tolerance ‘tol’ and a maximum number of iterations ‘maxiter.’

Next, we select the parameters the optimizer will drive by calling add_param and giving it the IndepVarComp unknowns that we have created. We also set high and low bounds for this problem. It is not required to set these (they will default to -1e99 and 1e99 respectively), but it is generally a good idea.

Finally, we add the objective. You can use any unknown in your model as the objective.

Once we have called setup on the model, we can specify the initial conditions for the design variables just like we did with unconnected params.

Since SLSQP is a gradient based optimizer, OpenMDAO will call the linearize method on the Paraboloid while calculating the total gradient of the objective with respect to the two design variables. This is done automatically.

Finally, we made a change to the print statement so that we can print the objective and the parameters. This time, we get the value by keying into the problem instance (‘top’) with the full variable path to the quantities we want to see. This is equivalent to what was shown in the first tutorial.

Putting this all together, when we run the model, we get output that looks like this (note, the optimizer may print some things before this, depending on settings):

...
Minimum of -27.333333 found at (6.666667, -7.333333)


## Optimization of the Paraboloid with a Constraint¶

Finally, let’s take this optimization problem and add a constraint to it. Our constraint takes the form of an inequality we want to satisfy: x - y >= 15.

First, we need to add one more import to the beginning of our model.

from openmdao.api import ExecComp


We’ll use an ExecComp to represent our constraint in the model. An ExecComp is a shortcut that lets us easily create a component that defines a simple expression for us.

top = Problem()

root = top.root = Group()

# Constraint Equation

root.connect('p1.x', 'p.x')
root.connect('p2.y', 'p.y')
root.connect('p.x', 'con.x')
root.connect('p.y', 'con.y')

top.driver = ScipyOptimizer()
top.driver.options['optimizer'] = 'SLSQP'

top.setup()
top.run()

print('\n')
print('Minimum of %f found at (%f, %f)' % (top['p.f_xy'], top['p.x'], top['p.y']))


Here, we added an ExecComp named ‘con’ to represent part of our constraint inequality. Our constraint is “x - y >= 15”, so we have created an ExecComp that will evaluate the expression “x - y” and place that result into the unknown ‘con.c’. To complete the definition of the constraint, we also need to connect our ‘con’ expression to ‘x’ and ‘y’ on the paraboloid.

Finally, we need to tell the driver to use the unknown “con.c” as a constraint using the add_constraint method. This method takes the name of the variable and an “upper” or “lower” bound. Here we give it a lower bound of 15, which completes the inequality constraint “x - y >= 15”.

OpenMDAO also supports the specification of double sided constraints, so if you wanted to constrain x-y to lie on a band between 15 and 16 which is “16 > x-y > 15”, you would just do the following:

top.driver.add_constraint('con.c', lower=15.0, upper=16.0)


So now, putting it all together, we can run the model and get this:

...
Minimum of -27.083333 found at (7.166667, -7.833333)


A new optimum is found because the original one was infeasible (i.e., that design point violated the constraint equation).

Tags