Design Space Exploration
In this demonstration we show how we can leverage the OptProblem class from the standard_evaluator library to perform optimization studies that can be version controlled.
The idea is that we will want to run multiple optimization problems that are driven by JSOn files that we can version control. In a future step we will also manage the assembly information itself in a JSON file that can be version controlled, but for this demonstration we went down an easier path.
We first load all required libraries and methods.
import typing
import openmdao.api as om
import openmdao.core.system as oms
import pandas as pd
from standard_evaluator import get_interface, show_structure, create_problem, get_state, set_state, load_state, save_state
from standard_evaluator import load_assembly, save_assembly, convert_om_var, get_opt_problem, set_opt_problem
from standard_evaluator import Variable, FloatVariable, IntVariable, ArrayVariable, OptProblem
matlabengine not available.
Creating the OpenMDAO assembly
In this example we use the Sellar Problem from the OpenMDAO documentation.
We first setup the assembly, set the optimization variables, constraints, and the objective. We can then run the optimization.
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'])
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.18339395172513
Iterations: 6
Function evaluations: 6
Gradient evaluations: 6
Optimization Complete
-----------------------------------
minimum found at
1.4707790339389824e-14
[1.97763888e+00 1.53488630e-14]
minumum objective
3.18339395172513
In the next step we get the description of the optimization problem and store that in an instance of the neutral OptProblem. We also capture the structure of the assembly in the neutral format.
We show the content of the original_opt variable,
original_opt = get_opt_problem(prob)
info = get_interface(prob.model)
original_opt
---------------------------------------------------------------------------
ValidationError Traceback (most recent call last)
Cell In[3], line 1
----> 1 original_opt = get_opt_problem(prob)
2 info = get_interface(prob.model)
3 original_opt
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/standard_evaluator/om_converter.py:607, in get_opt_problem(om_component)
605 responses = []
606 for info in local_dict['design_vars']:
--> 607 vars.append(convert_om_opt_vars(info[1]))
608 constraints = []
609 for info in local_dict['constraints']:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/standard_evaluator/om_converter.py:588, in convert_om_opt_vars(info)
585 input_dict['options'] = additional_dict
586 if info['size'] == 1:
587 # This is a float variable. Store it as such
--> 588 output = FloatVariable.model_validate(input_dict)
589 else:
590 # This is an array variable, and we need to capture its shape and
591 # additional information
592 print("Warning: More work is needed to support array variables")
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/pydantic/main.py:732, in BaseModel.model_validate(cls, obj, strict, extra, from_attributes, context, by_alias, by_name)
726 if by_alias is False and by_name is not True:
727 raise PydanticUserError(
728 'At least one of `by_alias` or `by_name` must be set to True.',
729 code='validate-by-alias-and-name-false',
730 )
--> 732 return cls.__pydantic_validator__.validate_python(
733 obj,
734 strict=strict,
735 extra=extra,
736 from_attributes=from_attributes,
737 context=context,
738 by_alias=by_alias,
739 by_name=by_name,
740 )
ValidationError: 2 validation errors for FloatVariable
bounds.0
Input should be a valid number [type=float_type, input_value=array([0.]), input_type=ndarray]
For further information visit https://errors.pydantic.dev/2.13/v/float_type
bounds.1
Input should be a valid number [type=float_type, input_value=array([10.]), input_type=ndarray]
For further information visit https://errors.pydantic.dev/2.13/v/float_type
Setting up a series of optimization studies
In the next step we take the original optimization problem and create a series of optimization problems where we change the upper bound of the first nonlinear constraint. This is a typical scenario in industry, we are trying to explore the design space and often the constraints are not fully defined.
We store all of the optimization problems in JSON files. This will allow us to fully describe the studies we are doing, and they can be version controlled.
In this example we are actually varying the bounds for both constraints, and not just one.
import numpy as np
# Define the lower and upper bounds
x_lower, x_upper = -2, 2
y_lower, y_upper = -22.5, -17
# Define the number of points along each axis
num_points_x = 5
num_points_y = 5
# Create evenly spaced values for x and y
x_values = np.linspace(x_lower, x_upper, num_points_x)
y_values = np.linspace(y_lower, y_upper, num_points_y)
# Create a meshgrid from the x and y values
X, Y = np.meshgrid(x_values, y_values)
# Combine the X and Y coordinates into a single array of points
grid_points = np.column_stack((X.ravel(), Y.ravel()))
# Add the point [0, 0] if it's not already included
#if not np.array_equal([0, 0], [x_lower, y_lower]):
# grid_points = np.vstack((grid_points, [0, 0]))
print(grid_points)
[[ -2. -22.5 ]
[ -1. -22.5 ]
[ 0. -22.5 ]
[ 1. -22.5 ]
[ 2. -22.5 ]
[ -2. -21.125]
[ -1. -21.125]
[ 0. -21.125]
[ 1. -21.125]
[ 2. -21.125]
[ -2. -19.75 ]
[ -1. -19.75 ]
[ 0. -19.75 ]
[ 1. -19.75 ]
[ 2. -19.75 ]
[ -2. -18.375]
[ -1. -18.375]
[ 0. -18.375]
[ 1. -18.375]
[ 2. -18.375]
[ -2. -17. ]
[ -1. -17. ]
[ 0. -17. ]
[ 1. -17. ]
[ 2. -17. ]]
import json
indx = 0
for upper1, upper2 in grid_points:
print(upper1, upper2)
new_opt = original_opt.model_copy()
new_opt.responses[0].bounds = (-np.inf, upper1)
new_opt.responses[1].bounds = (-np.inf, upper2)
optproblem_name = f'optproblem{indx:02}.json'
indx += 1
# Convert and write JSON object to file
with open(optproblem_name, "w") as outfile:
json.dump(new_opt.model_dump(exclude_unset=True), outfile, indent=2, skipkeys=True)
-2.0 -22.5
-1.0 -22.5
0.0 -22.5
1.0 -22.5
2.0 -22.5
-2.0 -21.125
-1.0 -21.125
0.0 -21.125
1.0 -21.125
2.0 -21.125
-2.0 -19.75
-1.0 -19.75
0.0 -19.75
1.0 -19.75
2.0 -19.75
-2.0 -18.375
-1.0 -18.375
0.0 -18.375
1.0 -18.375
2.0 -18.375
-2.0 -17.0
-1.0 -17.0
0.0 -17.0
1.0 -17.0
2.0 -17.0
Running a series of optimization problems
In the next step we actually run the different optimization problems. To do this we first create a pattern for all the optimization files, and get the names of all of them.
Then we run a loop through all the file names, first loading the optimization problem formulation into a dictionary, convert it into OptProblem’s, create a new instance of the Sellar problem, set the optimization problem, the optimzition parameters, and run the optimization. We store the results from the optimization in a Pandas DataFrame together with the name of the file that defined the specific optimization problem.
After running all the optimizations we combine the results into a single DataFrame, and show the summary.
solutions = []
import glob
# Specify the pattern for the JSON files
pattern = 'optproblem*.json'
# Use glob to find all files matching the pattern
json_files = glob.glob(pattern)
# Loop through the list of files and load each JSON file
data_list = []
for file in json_files:
with open(file, 'r') as f:
data = json.load(f) # Load the JSON data
new_opt = OptProblem.model_validate(data)
new_prob = om.Problem()
new_prob.model = SellarMDA()
# Setting the optimization problem
set_opt_problem(new_prob, new_opt)
# Setting the optimizer options.
new_prob.driver = om.ScipyOptimizeDriver()
new_prob.driver.options['optimizer'] = 'SLSQP'
# prob.driver.options['maxiter'] = 100
new_prob.driver.options['tol'] = 1e-8
# Ask OpenMDAO to finite-difference across the model to compute the gradients for the optimizer
new_prob.model.approx_totals()
new_prob.set_solver_print(level=0)
# Run the optimization
new_prob.run_driver()
# Check if the optimization was successful
if not new_prob.driver.fail:
print("Optimization was successful!")
# Get the state of the assembly, which will be the optimization solution
state = get_state(new_prob, info)
# We will store the solution in a dictionary, including the name of the file.
state = {key: value.flatten()[0] if len(value) == 1 else value for key, value in state.items()}
state['name'] = file
# Convert the dictionary to a Pandas Series, and store as a DataFrame
solutions.append(pd.Series(state).to_frame().T)
else:
print("Optimization failed.")
# Combine all the solutions into a single DataFrame
df = pd.concat(solutions)
# Show the result DataFrame
df
Positive directional derivative for linesearch (Exit mode 8)
Current function value: 36.7744422015937
Iterations: 16
Function evaluations: 62
Gradient evaluations: 12
Optimization FAILED.
Positive directional derivative for linesearch
-----------------------------------
Optimization failed.
Positive directional derivative for linesearch (Exit mode 8)
Current function value: 25.106465928943965
Iterations: 11
Function evaluations: 17
Gradient evaluations: 7
Optimization FAILED.
Positive directional derivative for linesearch
-----------------------------------
Optimization failed.
NL: NLBGSSolver 'NL: NLBGS' on system 'cycle' failed to converge in 10 iterations.
Inequality constraints incompatible (Exit mode 4)
Current function value: 15.314007967433286
Iterations: 13
Function evaluations: 58
Gradient evaluations: 12
Optimization FAILED.
Inequality constraints incompatible
-----------------------------------
Optimization failed.
Optimization terminated successfully (Exit mode 0)
Current function value: 8.316848498193089
Iterations: 9
Function evaluations: 8
Gradient evaluations: 8
Optimization Complete
-----------------------------------
Optimization was successful!
NL: NLBGSSolver 'NL: NLBGS' on system 'cycle' failed to converge in 10 iterations.
NL: NLBGSSolver 'NL: NLBGS' on system 'cycle' failed to converge in 10 iterations.
Optimization terminated successfully (Exit mode 0)
Current function value: 2.8815345618132655
Iterations: 6
Function evaluations: 7
Gradient evaluations: 6
Optimization Complete
-----------------------------------
Optimization was successful!
Optimization terminated successfully (Exit mode 0)
Current function value: 32.15279509345695
Iterations: 6
Function evaluations: 6
Gradient evaluations: 6
Optimization Complete
-----------------------------------
Optimization was successful!
Optimization terminated successfully (Exit mode 0)
Current function value: 20.258749369476174
Iterations: 6
Function evaluations: 6
Gradient evaluations: 6
Optimization Complete
-----------------------------------
Optimization was successful!
Optimization terminated successfully (Exit mode 0)
Current function value: 11.27091613374096
Iterations: 6
Function evaluations: 6
Gradient evaluations: 6
Optimization Complete
-----------------------------------
Optimization was successful!
Optimization terminated successfully (Exit mode 0)
Current function value: 5.389808017048902
Iterations: 6
Function evaluations: 7
Gradient evaluations: 6
Optimization Complete
-----------------------------------
Optimization was successful!
NL: NLBGSSolver 'NL: NLBGS' on system 'cycle' failed to converge in 10 iterations.
Optimization terminated successfully (Exit mode 0)
Current function value: 2.615642164038142
Iterations: 8
Function evaluations: 8
Gradient evaluations: 8
Optimization Complete
-----------------------------------
Optimization was successful!
Optimization terminated successfully (Exit mode 0)
Current function value: 9.566593149581722
Iterations: 7
Function evaluations: 7
Gradient evaluations: 7
Optimization Complete
-----------------------------------
Optimization was successful!
Optimization terminated successfully (Exit mode 0)
Current function value: 4.189681516258113
Iterations: 7
Function evaluations: 7
Gradient evaluations: 7
Optimization Complete
-----------------------------------
Optimization was successful!
Optimization terminated successfully (Exit mode 0)
Current function value: 3.183393952145962
Iterations: 7
Function evaluations: 7
Gradient evaluations: 7
Optimization Complete
-----------------------------------
Optimization was successful!
Optimization terminated successfully (Exit mode 0)
Current function value: 2.2033093084417064
Iterations: 8
Function evaluations: 8
Gradient evaluations: 8
Optimization Complete
-----------------------------------
Optimization was successful!
NL: NLBGSSolver 'NL: NLBGS' on system 'cycle' failed to converge in 10 iterations.
NL: NLBGSSolver 'NL: NLBGS' on system 'cycle' failed to converge in 10 iterations.
Optimization terminated successfully (Exit mode 0)
Current function value: 1.2549821858007038
Iterations: 9
Function evaluations: 9
Gradient evaluations: 9
Optimization Complete
-----------------------------------
Optimization was successful!
Optimization terminated successfully (Exit mode 0)
Current function value: 5.168711366061382
Iterations: 6
Function evaluations: 6
Gradient evaluations: 6
Optimization Complete
-----------------------------------
Optimization was successful!
Optimization terminated successfully (Exit mode 0)
Current function value: 4.173853524365584
Iterations: 7
Function evaluations: 7
Gradient evaluations: 7
Optimization Complete
-----------------------------------
Optimization was successful!
Optimization terminated successfully (Exit mode 0)
Current function value: 3.1833939516764156
Iterations: 6
Function evaluations: 6
Gradient evaluations: 6
Optimization Complete
-----------------------------------
Optimization was successful!
Optimization terminated successfully (Exit mode 0)
Current function value: 2.2033093084548656
Iterations: 6
Function evaluations: 7
Gradient evaluations: 6
Optimization Complete
-----------------------------------
Optimization was successful!
NL: NLBGSSolver 'NL: NLBGS' on system 'cycle' failed to converge in 10 iterations.
NL: NLBGSSolver 'NL: NLBGS' on system 'cycle' failed to converge in 10 iterations.
Optimization terminated successfully (Exit mode 0)
Current function value: 1.2549821858072014
Iterations: 7
Function evaluations: 7
Gradient evaluations: 7
Optimization Complete
-----------------------------------
Optimization was successful!
Optimization terminated successfully (Exit mode 0)
Current function value: 5.168711365989337
Iterations: 6
Function evaluations: 6
Gradient evaluations: 6
Optimization Complete
-----------------------------------
Optimization was successful!
Optimization terminated successfully (Exit mode 0)
Current function value: 4.173853524378517
Iterations: 6
Function evaluations: 6
Gradient evaluations: 6
Optimization Complete
-----------------------------------
Optimization was successful!
Optimization terminated successfully (Exit mode 0)
Current function value: 3.183393951729169
Iterations: 6
Function evaluations: 6
Gradient evaluations: 6
Optimization Complete
-----------------------------------
Optimization was successful!
Optimization terminated successfully (Exit mode 0)
Current function value: 2.2033093084479907
Iterations: 6
Function evaluations: 7
Gradient evaluations: 6
Optimization Complete
-----------------------------------
Optimization was successful!
NL: NLBGSSolver 'NL: NLBGS' on system 'cycle' failed to converge in 10 iterations.
NL: NLBGSSolver 'NL: NLBGS' on system 'cycle' failed to converge in 10 iterations.
Optimization terminated successfully (Exit mode 0)
Current function value: 1.2549821858054697
Iterations: 7
Function evaluations: 7
Gradient evaluations: 7
Optimization Complete
-----------------------------------
Optimization was successful!
Optimization terminated successfully (Exit mode 0)
Current function value: 3.183393951729169
Iterations: 6
Function evaluations: 6
Gradient evaluations: 6
Optimization Complete
-----------------------------------
Optimization was successful!
| y1 | y2 | z | x | con1 | con2 | obj | name | |
|---|---|---|---|---|---|---|---|---|
| 0 | 2.16 | 1.5 | [0.0, 0.03030615433010159] | 2.429694 | 1.0 | -22.5 | 8.316848 | optproblem03.json |
| 0 | 1.16 | 1.5 | [1.056894812204979e-13, 0.422967038573036] | 1.037033 | 2.0 | -22.5 | 2.881535 | optproblem04.json |
| 0 | 5.16 | 2.875 | [0.0, 0.6034366616784598] | 5.131563 | -2.0 | -21.125 | 32.152795 | optproblem05.json |
| 0 | 4.16 | 2.875 | [6.582959922515623e-15, 0.835392194558752] | 3.899608 | -1.0 | -21.125 | 20.258749 | optproblem06.json |
| 0 | 3.16 | 2.875 | [0.0, 1.097361116516697] | 2.637639 | 0.0 | -21.125 | 11.270916 | optproblem07.json |
| 0 | 2.16 | 2.875 | [0.0, 1.4053061543301015] | 1.329694 | 1.0 | -21.125 | 5.389808 | optproblem08.json |
| 0 | 1.16 | 2.109774 | [0.0, 1.0327406260261245] | 0.549214 | 2.0 | -21.890226 | 2.615642 | optproblem09.json |
| 0 | 5.16 | 4.25 | [1.978436661679349, 1.34947960803124e-12] | 2.095788 | -2.0 | -19.75 | 9.566593 | optproblem10.json |
| 0 | 4.16 | 4.25 | [2.210392194518491, 4.795678246750001e-13] | 0.124166 | -1.0 | -19.75 | 4.189682 | optproblem11.json |
| 0 | 3.16 | 3.755278 | [1.977638883605749, 7.444231462403437e-14] | 0.0 | -0.0 | -20.244722 | 3.183394 | optproblem12.json |
| 0 | 2.16 | 3.139388 | [1.6696938456669281, 3.8634046702380526e-14] | 0.0 | 1.0 | -20.860612 | 2.203309 | optproblem13.json |
| 0 | 1.16 | 2.354066 | [1.2770329614250442, 1.4182406177508744e-14] | 0.0 | 2.0 | -21.645934 | 1.254982 | optproblem14.json |
| 0 | 5.16 | 4.743127 | [2.471563338346374, 3.851586669543894e-12] | 0.0 | -2.0 | -19.256873 | 5.168711 | optproblem15.json |
| 0 | 4.16 | 4.279216 | [2.2396078054336765, 1.5650485625672098e-13] | 0.0 | -1.0 | -19.720784 | 4.173854 | optproblem16.json |
| 0 | 3.16 | 3.755278 | [1.9776388834729854, 1.8925849432621544e-15] | 0.0 | -0.0 | -20.244722 | 3.183394 | optproblem17.json |
| 0 | 2.16 | 3.139388 | [1.6696938456715893, 8.66652570560973e-15] | 0.0 | 1.0 | -20.860612 | 2.203309 | optproblem18.json |
| 0 | 1.16 | 2.354066 | [1.2770329614291231, 1.951237785459953e-15] | 0.0 | 2.0 | -21.645934 | 1.254982 | optproblem19.json |
| 0 | 5.16 | 4.743127 | [2.4715633383321443, 0.0] | 0.0 | -2.0 | -19.256873 | 5.168711 | optproblem20.json |
| 0 | 4.16 | 4.279216 | [2.239607805437483, 3.363869713064246e-14] | 0.0 | -1.0 | -19.720784 | 4.173854 | optproblem21.json |
| 0 | 3.16 | 3.755278 | [1.977638883487764, 8.830566052859473e-15] | 0.0 | -0.0 | -20.244722 | 3.183394 | optproblem22.json |
| 0 | 2.16 | 3.139388 | [1.6696938456710986, 1.4729587173841369e-15] | 0.0 | 1.0 | -20.860612 | 2.203309 | optproblem23.json |
| 0 | 1.16 | 2.354066 | [1.2770329614285358, 3.876178038906271e-15] | 0.0 | 2.0 | -21.645934 | 1.254982 | optproblem24.json |
| 0 | 3.16 | 3.755278 | [1.977638883487764, 8.830566052859473e-15] | 0.0 | -0.0 | -20.244722 | 3.183394 | optproblem25.json |
Since in this example we only varied one constraint we can create a plot showing the impact of the change of the constraint on the value of the objective function, which is what we show next to finish this demonstration.
import matplotlib.pyplot as plt
# Ensure that plots are displayed inline in Jupyter Notebook
%matplotlib inline
# Create a 3D plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
# Plot the data
ax.scatter(df['con1'], df['con2'], df['obj'], c='b', marker='o') # 'o' for dots
# Set labels
ax.set_xlabel('con1')
ax.set_ylabel('con2')
ax.set_zlabel('obj')
# Show the plot
plt.show()
