Advanced Recording Example

Below we demonstrate a more advanced example of case recording including the four different objects that a recorder can be attached to. We will then show how to extract various data from the model and finally, relate that back to an XDSM for illustrative purposes.

import numpy as np
from openmdao.test_suite.components.sellar_feature import SellarMDAWithUnits

# build the model
prob = om.Problem(model=SellarMDAWithUnits())

model = prob.model
model.add_design_var('z', lower=np.array([-10.0, 0.0]),
                          upper=np.array([10.0, 10.0]))
model.add_design_var('x', lower=0.0, upper=10.0)
model.add_objective('obj')
model.add_constraint('con1', upper=0.0)
model.add_constraint('con2', upper=0.0)

# setup the optimization
prob.driver = om.ScipyOptimizeDriver(optimizer='SLSQP', tol=1e-9, disp=False)

# Here we show how to attach recorders to each of the four objects:
#   problem, driver, solver, and system

# Create a recorder
recorder = om.SqliteRecorder('cases.sql')

# Attach recorder to the problem
prob.add_recorder(recorder)

# Attach recorder to the driver
prob.driver.add_recorder(recorder)

# To attach a recorder to a subsystem or solver, you need to call `setup`
# first so that the model hierarchy has been generated
prob.setup()

# Attach recorder to a subsystem
model.obj_cmp.add_recorder(recorder)

# Attach recorder to a solver
model.cycle.nonlinear_solver.add_recorder(recorder)

prob.set_solver_print(0)
prob.run_driver()
prob.record("final_state")
prob.cleanup()

The following XDSM diagram shows the SellarMDA component equations and their inputs and outputs. Through the different recorders we can access the different parts of the model. We’ll take you through an example of each object and relate it back to the diagram.

sellar_xdsm

Driver

If we want to view the convergence of the model, the best place to find that by looking at the cases recorded by the driver. By default, a recorder attached to a driver will record the design variables, constraints and objectives, so we will print them for the model at the end of the optimization. We’ll use the helper methods like get_objectives, get_design_vars, get_constraints to return the info we’re seeking.

import openmdao.api as om

# Instantiate your CaseReader
cr = om.CaseReader("cases.sql")

# List driver cases (do not recurse to system/solver cases)
driver_cases = cr.list_cases('driver', recurse=False);
last_case = cr.get_case(driver_cases[-1])

objectives = last_case.get_objectives()
design_vars = last_case.get_design_vars()
constraints = last_case.get_constraints()

print(objectives['obj'])
print(design_vars['x'])
print(design_vars['z'])
print(constraints['con1'])
print(constraints['con2'])
driver
rank0:ScipyOptimize_SLSQP|0
rank0:ScipyOptimize_SLSQP|1
rank0:ScipyOptimize_SLSQP|2
rank0:ScipyOptimize_SLSQP|3
rank0:ScipyOptimize_SLSQP|4
rank0:ScipyOptimize_SLSQP|5
rank0:ScipyOptimize_SLSQP|6
rank0:ScipyOptimize_SLSQP|7
rank0:ScipyOptimize_SLSQP|8
rank0:ScipyOptimize_SLSQP|9
rank0:ScipyOptimize_SLSQP|10
rank0:ScipyOptimize_SLSQP|11
rank0:ScipyOptimize_SLSQP|12
rank0:ScipyOptimize_SLSQP|13
[3.18339395]
[0.]
[1.97763888e+00 1.25035459e-15]
[-1.68550507e-10]
[-20.24472223]

Plotting Design Variables

When inspecting or debugging a model, it can be helpful to visualize the path of the design variables to their final values. To do this, we can list the cases of the driver and plot the data with respect to the iteration number.

import matplotlib.pyplot as plt
import numpy as np

# Instantiate your CaseReader
cr = om.CaseReader("cases.sql")

# List driver cases (do not recurse to system/solver cases, suppress display)
driver_cases = cr.list_cases('driver', recurse=False, out_stream=None)

# Plot the path the design variables took to convergence
# Note that there are two lines in the right plot because "Z"
# contains two variables that are being optimized
dv_x_values = []
dv_z_values = []
for case_id in driver_cases:
    case = cr.get_case(case_id)
    design_vars = case.get_design_vars()
    dv_x_values.append(design_vars['x'])
    dv_z_values.append(design_vars['z'])

fig, (ax1, ax2) = plt.subplots(1, 2)
fig.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.5, hspace=None)
ax1.plot(np.arange(len(dv_x_values)), np.array(dv_x_values))

ax1.set(xlabel='Iterations', ylabel='Design Var: X', title='Optimization History')
ax1.grid()

ax2.plot(np.arange(len(dv_z_values)), np.array(dv_z_values))

ax2.set(xlabel='Iterations', ylabel='Design Var: Z', title='Optimization History')
ax2.grid()

plt.show()
../../_images/advanced_case_recording_8_0.png

Problem

A Problem recorder is best if you want to record an arbitrary case before or after a running the model. In our case, we have placed our point at the end of the model.

sorted(problem_vars['outputs'])
['con1', 'con2', 'obj', 'x', 'y1', 'y2', 'z']

System

Suppose we want to know the value of y1 going into the objective function. For this we will use the system recorder. Using the list_cases method, we’ll get the list of cases that were recorded by the objective component root.obj_comp. Then we use get_case to determine the input names using the keys of the first case’s inputs dictionary. Since we want to find the value of y1 going into the objective function for each execution, we use the get_cases function to iterate through those 14 cases and get the value in each.

import numpy as np

# Instantiate your CaseReader
cr = om.CaseReader("cases.sql")

system_cases = cr.list_cases('root.obj_cmp', recurse=False)

# Get the keys of all the inputs to the objective function
case = cr.get_case(system_cases[0])

cr.get_cases('root.obj_cmp');
root.obj_cmp
rank0:ScipyOptimize_SLSQP|0|root._solve_nonlinear|0|NLRunOnce|0|obj_cmp._solve_nonlinear|0
rank0:ScipyOptimize_SLSQP|1|root._solve_nonlinear|1|NLRunOnce|0|obj_cmp._solve_nonlinear|1
rank0:ScipyOptimize_SLSQP|2|root._solve_nonlinear|2|NLRunOnce|0|obj_cmp._solve_nonlinear|2
rank0:ScipyOptimize_SLSQP|3|root._solve_nonlinear|3|NLRunOnce|0|obj_cmp._solve_nonlinear|3
rank0:ScipyOptimize_SLSQP|4|root._solve_nonlinear|4|NLRunOnce|0|obj_cmp._solve_nonlinear|4
rank0:ScipyOptimize_SLSQP|5|root._solve_nonlinear|5|NLRunOnce|0|obj_cmp._solve_nonlinear|5
rank0:ScipyOptimize_SLSQP|6|root._solve_nonlinear|6|NLRunOnce|0|obj_cmp._solve_nonlinear|6
rank0:ScipyOptimize_SLSQP|7|root._solve_nonlinear|7|NLRunOnce|0|obj_cmp._solve_nonlinear|7
rank0:ScipyOptimize_SLSQP|8|root._solve_nonlinear|8|NLRunOnce|0|obj_cmp._solve_nonlinear|8
rank0:ScipyOptimize_SLSQP|9|root._solve_nonlinear|9|NLRunOnce|0|obj_cmp._solve_nonlinear|9
rank0:ScipyOptimize_SLSQP|10|root._solve_nonlinear|10|NLRunOnce|0|obj_cmp._solve_nonlinear|10
rank0:ScipyOptimize_SLSQP|11|root._solve_nonlinear|11|NLRunOnce|0|obj_cmp._solve_nonlinear|11
rank0:ScipyOptimize_SLSQP|12|root._solve_nonlinear|12|NLRunOnce|0|obj_cmp._solve_nonlinear|12
rank0:ScipyOptimize_SLSQP|13|root._solve_nonlinear|13|NLRunOnce|0|obj_cmp._solve_nonlinear|13

Solver

Similar to the system recorder, we can query the cases recorded by the solver. You can also access the values of inputs to the equation with the solver but in this case we’ll focus on the coupling variables y1 and y2.

We’ll use list_cases to get the cases recorded for 'root.cycle.nonlinear_solver', find out how many cases there are and plot a subset of the solver iterations (covering multiple driver iterations) to see the final convergence.

import numpy as np
import matplotlib.pyplot as plt
import openmdao.api as om

# Instantiate your CaseReader
cr = om.CaseReader("cases.sql")

# Get list of cases, without displaying them
solver_cases = cr.list_cases('root.cycle.nonlinear_solver', out_stream=None)

# Plot the convergence of the two coupling variables (last 35 iterations)
y1_history = []
y2_history = []

for case_id in solver_cases[-35:]:
    case = cr.get_case(case_id)
    y1_history.append(case['y1'])
    y2_history.append(case['y2'])

iterations = np.arange(-len(y1_history), 0, 1)

fig, (ax1, ax2) = plt.subplots(2, 1)

ax1.plot(iterations, np.array(y1_history))
ax1.set(ylabel='Coupling Output: y1', title='Solver History')
ax1.grid()

ax2.plot(iterations, np.array(y2_history))
ax2.set(ylabel='Coupling Parameter: y2', xlabel='Iterations')
ax2.grid()

plt.show()

# Get the final values
case = cr.get_case(solver_cases[-1])
print(case['y1'])
print(case['y2'])

print(f"Number of cases: {len(solver_cases)}")
../../_images/advanced_case_recording_17_0.png
[3.16]
[3.75527777]
Number of cases: 76