5. Optimization
The notebook will walk through a simple example of performing optimization of a CSDL model using modopt.
Begin by initializing a csdl Recorder object.
import warnings
warnings.filterwarnings('ignore')
import csdl_alpha as csdl
import numpy as np
recorder = csdl.Recorder(inline=True)
recorder.start()
In this example, we will use the toy sellar optimization problem as exactly defined in OpenMDAO’s documentation. The problem is a simple system with coupled disciplines. See their documentation page for more details.
We define the problem here using namespaces, nonlinear solvers, as well as special methods to specify optimization variables: set_as_design_variable, set_as_constraint and set_as_objective.
z1 = csdl.Variable(name = 'z1', value = 5)
z2 = csdl.Variable(name = 'z2', value = 2)
x = csdl.Variable(name = 'x', value = 1.0)
y2 = csdl.ImplicitVariable(name = 'y2', value = 1.0)
z1.set_as_design_variable(lower=-10.0, upper=10.0) # design variable
z2.set_as_design_variable(lower= 0.0, upper=10.0) # design variable
x.set_as_design_variable(lower= 0.0, upper=10.0) # design variable
# Define each "component" from the example
with csdl.namespace('Discipline 1'):
y1 = z1**2 + z2 + x - 0.2*y2
y1.add_name('y1')
with csdl.namespace('Discipline 2'):
residual = csdl.sqrt(y1) + z1 + z2 - y2
residual.add_name('residual')
with csdl.namespace('Objective'):
f = x**2 + z2 + y1 + csdl.exp(-y2)
f.add_name('f')
f.set_as_objective() # objective
with csdl.namespace('Constraint 1'):
g1 = 3.16 - y1
g1.add_name('g1')
g1.set_as_constraint(upper=0.0) # constraint
with csdl.namespace('Constraint 2'):
g2 = y2 - 24.0
g2.add_name('g2')
g2.set_as_constraint(upper=0.0) # constraint
# Specifiy coupling
with csdl.namespace('Couple'):
solver = csdl.nonlinear_solvers.Newton()
solver.add_state(y2, residual, tolerance=1e-8)
solver.run()
recorder.stop()
nonlinear solver: newton_nlsolver converged in 2 iterations.
Here we set three design variables, two constraints and an objective. Because we instantiated the Recorder with inline = True, we can see the values of the initial evaluation:
print('intital run values:')
print('y1: ', y1.value)
print('y2: ', y2.value)
print('f: ', f.value)
print('g1: ', g1.value)
print('g2: ', g2.value)
print('z1: ', z1.value)
print('z2: ', z2.value)
intital run values:
y1: [25.58830237]
y2: [12.05848815]
f: [28.58830816]
g1: [-22.42830237]
g2: [-11.94151185]
z1: [5.]
z2: [2.]
CSDL makes it easy to evaluate the model by providing the Simulator class. As the name suggests, the Simulator class is an interface for evaluating your model and any derivatives after the model is defined.
We currently provide two experimental simulators: PySimulator and JaxSimulator which computes your model using numpy and jax respectively. For this example, we will use the PySimulator for easier use. See the ‘JAX’ tutorial to learn more about running your model efficiently with JaxSimulator. Creating a PySimulator is as simple as instantiating it with your Recorder object.
sim = csdl.experimental.PySimulator(recorder)
To interface our simulator with modopt, we use modopt’s CSDLAlphaProblem. We can then solve the optimization using modopt’s API.
from modopt import CSDLAlphaProblem
from modopt import SLSQP
# Instantiate your problem using the csdl Simulator object and name your problem
prob = CSDLAlphaProblem(problem_name='sellar',simulator=sim)
optimizer = SLSQP(prob, ftol=1e-9, maxiter=20)
# Check first derivatives at the initial guess, if needed
# optimizer.check_first_derivatives(prob.x0)
# Solve your optimization problem
optimizer.solve()
nonlinear solver: newton_nlsolver converged in 2 iterations.
nonlinear solver: newton_nlsolver converged in 2 iterations.
Deleting self.pFpx ...
Deleting self.pCpx, pCpx_dict ...
Computing constraints >>>>>>>>>>
---------Computed constraints---------
Computing objective >>>>>>>>>>
---------Computed objective---------
Computing gradient >>>>>>>>>>
---------Computed gradient---------
Computing constraints >>>>>>>>>>
---------Computed constraints---------
Computing Jacobian >>>>>>>>>>
---------Computed Jacobian---------
Computing objective >>>>>>>>>>
nonlinear solver: newton_nlsolver converged in 2 iterations.
---------Computed objective---------
Computing constraints >>>>>>>>>>
---------Computed constraints---------
Computing gradient >>>>>>>>>>
nonlinear solver: newton_nlsolver converged in 2 iterations.
---------Computed gradient---------
Computing Jacobian >>>>>>>>>>
---------Computed Jacobian---------
Computing objective >>>>>>>>>>
nonlinear solver: newton_nlsolver converged in 3 iterations.
---------Computed objective---------
Computing constraints >>>>>>>>>>
---------Computed constraints---------
Computing gradient >>>>>>>>>>
nonlinear solver: newton_nlsolver converged in 3 iterations.
---------Computed gradient---------
Computing Jacobian >>>>>>>>>>
---------Computed Jacobian---------
Computing objective >>>>>>>>>>
nonlinear solver: newton_nlsolver converged in 3 iterations.
---------Computed objective---------
Computing constraints >>>>>>>>>>
---------Computed constraints---------
Computing gradient >>>>>>>>>>
nonlinear solver: newton_nlsolver converged in 3 iterations.
---------Computed gradient---------
Computing Jacobian >>>>>>>>>>
---------Computed Jacobian---------
Computing objective >>>>>>>>>>
nonlinear solver: newton_nlsolver converged in 3 iterations.
---------Computed objective---------
Computing constraints >>>>>>>>>>
---------Computed constraints---------
Computing gradient >>>>>>>>>>
nonlinear solver: newton_nlsolver converged in 3 iterations.
---------Computed gradient---------
Computing Jacobian >>>>>>>>>>
---------Computed Jacobian---------
Computing objective >>>>>>>>>>
nonlinear solver: newton_nlsolver converged in 3 iterations.
---------Computed objective---------
Computing constraints >>>>>>>>>>
---------Computed constraints---------
Computing gradient >>>>>>>>>>
nonlinear solver: newton_nlsolver converged in 3 iterations.
---------Computed gradient---------
Computing Jacobian >>>>>>>>>>
---------Computed Jacobian---------
Finally, lets print the optimized values and verify our results match OpenMDAO’s results.
optimizer.print_results()
print('optimized values:')
print('y1: ', y1.value)
print('y2: ', y2.value)
print('f: ', f.value)
print('g1: ', g1.value)
print('g2: ', g2.value)
print('z1: ', z1.value)
print('z2: ', z2.value)
assert np.isclose(y1.value, 3.15999999)
assert np.isclose(y2.value, 3.75527776)
assert np.isclose(f.value, 3.18339394)
assert np.isclose(g1.value, -8.95310492e-11)
assert np.isclose(g2.value, -20.24472224)
assert np.isclose(z1.value, 1.97763888)
assert np.isclose(z2.value, 0)
==============
Scipy summary:
==============
Problem : sellar
Solver : scipy_slsqp
Success : True
Message : Optimization terminated successfully
Objective : 3.183393951728967
Gradient norm : 4.022879450485861
Total time : 0.01696324348449707
Major iterations : 6
Total function evals : 6
Total gradient evals : 6
==================================================
optimized values:
y1: [3.16]
y2: [3.75527777]
f: [3.18339395]
g1: [-8.95310492e-11]
g2: [-20.24472223]
z1: [1.97763888]
z2: [0.]
We often want to run the model stand-alone without any optimization. Simulator provides this functionality by the run and compute_totals methods. Use the getters and setters to see and change the values of variables.
print('current f: ', sim[f])
print('current dfdz1: ', sim.compute_totals(f,z1)[f,z1], '\n')
sim[z1] += 1.0
sim.run()
dfdz1 = sim.compute_totals(f,z1)[f,z1]
print('\nnew f: ', sim[f])
print('current dfdz1: ', dfdz1)
current f: [3.18339395]
nonlinear solver: newton_nlsolver converged in 3 iterations.
current dfdz1: [[3.50848986]]
nonlinear solver: newton_nlsolver converged in 2 iterations.
nonlinear solver: newton_nlsolver converged in 2 iterations.
new f: [7.71844379]
current dfdz1: [[5.54894571]]
Use check_totals to verify analytical derivatives using finite difference.
dfdz1 = sim.check_totals()
nonlinear solver: newton_nlsolver converged in 2 iterations.
nonlinear solver: newton_nlsolver converged in 2 iterations.
nonlinear solver: newton_nlsolver converged in 2 iterations.
nonlinear solver: newton_nlsolver converged in 2 iterations.
nonlinear solver: newton_nlsolver converged in 2 iterations.
nonlinear solver: newton_nlsolver converged in 2 iterations.
Derivative Verification Results
-------------------------------
ofs (3) wrts (3) norm fd norm error rel error tags
------------------------------------------------------------------------------------------------------------------------------
f z1 5.548945712851889 5.548946721667392 1.0088155022458523e-06 1.8180306152637114e-07 obj(1,),dv(1,),
f z2 1.7685933165868333 1.7685933180189295 1.4320962193892228e-09 8.097374364126682e-10 obj(1,),dv(1,),
f x 0.9646991795541301 0.9647001810719757 1.0015178456557194e-06 1.0381648778616708e-06 obj(1,),dv(1,),
g1 z1 5.555277766964987 5.5552787694423955 1.0024774086403454e-06 1.8045492409032929e-07 con(1,),dv(1,),
g1 z2 0.7721994304206554 0.7721994297327228 6.87932599774399e-10 8.908742654893042e-10 con(1,),dv(1,),
g1 x 0.9652492880258193 0.9652492893863496 1.360530244021163e-09 1.4095117799942754e-09 con(1,),dv(1,),
g2 z1 2.0000000000580482 1.9999999913977717 8.66027649593093e-09 4.330138266589884e-09 con(1,),dv(1,),
g2 z2 1.1390028478967227 1.1390028475943836 3.023390426903916e-10 2.654418672692021e-10 con(1,),dv(1,),
g2 x 0.17375355987090357 0.17375354843807145 1.1432832125457537e-08 6.57991288709271e-08 con(1,),dv(1,),