Divya Manglam
Reputation: 491
How to implement SAND architecture using OpenMDAO
I am trying to implement Simultaneous ANalysis and Design (SAND) architecture on OpenMDAO with Sellar Problem. I thought of following way to do so -
class SellarDis1(Component):
def __init__(self):
super(SellarDis1, self).__init__()
self.add_param('z', val=np.zeros(2))
self.add_param('x', val=0.0)
self.add_param('y2', val=1.0)
self.add_output('y1', val=1.0)
def solve_nonlinear(self, params, unknowns, resids):
pass
def apply_nonlinear(self, params, unknowns, resids):
z1 = params['z'][0]
z2 = params['z'][1]
x1 = params['x']
y2 = params['y2']
y1 = unknowns['y1']
resids['y1'] = z1**2 + z2 + x1 - 0.2*y2 - y1
def linearize(self, params, unknowns, resids):
J = {}
J['y1','y1'] = 1.0
J['y1','y2'] = -0.2
J['y1','z'] = np.array([[2*params['z'][0], 1.0]])
J['y1','x'] = 1.0
return J
class SellarDis2(Component):
def __init__(self):
super(SellarDis2, self).__init__()
self.add_param('z', val=np.zeros(2))
self.add_param('y1', val=1.0)
self.add_output('y2', val=1.0)
def solve_nonlinear(self, params, unknowns, resids):
pass
def apply_nonlinear(self, params, unknowns, resids):
z1 = params['z'][0]
z2 = params['z'][1]
y1 = params['y1']
y1 = abs(y1)
y2 = unknowns['y2']
resids['y2'] = y1**.5 + z1 + z2 - y2
def linearize(self, params, unknowns, resids):
J = {}
J['y2', 'y2'] = 1.0
J['y2', 'y1'] = 0.5*params['y1']**-0.5
J['y2', 'z'] = np.array([[1.0, 1.0]])
return J
class SellarSAND(Group):
def __init__(self):
super(SellarSAND, self).__init__()
self.add('px', IndepVarComp('x', 1.0), promotes=['*'])
self.add('pz', IndepVarComp('z', np.array([5.0,2.0])), promotes=['*'])
self.add('d1', SellarDis1(), promotes=['*'])
self.add('d2', SellarDis2(), promotes=['*'])
self.add('obj_cmp', ExecComp('obj = x**2 + z[1] + y1 + exp(-y2)',
z=np.array([0.0, 0.0]), x=0.0, y1=0.0, y2=0.0),
promotes=['*'])
self.add('con_cmp1', ExecComp('con1 = 3.16 - y1'), promotes=['*'])
self.add('con_cmp2', ExecComp('con2 = y2 - 24.0'), promotes=['*'])
self.nl_solver = NLGaussSeidel()
self.nl_solver.options['atol'] = 1.0e-12
self.ln_solver = ScipyGMRES()
top = Problem()
top.root = SellarSAND()
top.driver = ScipyOptimizer()
top.driver.options['optimizer'] = 'SLSQP'
top.driver.options['tol'] = 1.0e-12
top.driver.add_desvar('z', lower=np.array([-10.0, 0.0]),upper=np.array([10.0, 10.0]))
top.driver.add_desvar('x', lower=0.0, upper=10.0)
top.driver.add_objective('obj')
top.driver.add_constraint('con1', upper=0.0)
top.driver.add_constraint('con2', upper=0.0)
top.setup()
tt = time.time()
top.run()
print("\n")
print( "Minimum found at (z1,z2,x) = (%f, %f, %f)" % (top['z'][0], \
top['z'][1], \
top['x']))
print("Coupling vars: %f, %f" % (top['y1'], top['y2']))
print("Minimum objective: ", top['obj'])
print("Elapsed time: ", 1000*(time.time()-tt), "milliseconds")
But on running the system I am getting following incorrect result -
Iteration limit exceeded (Exit mode 9)
Current function value: [ 1.36787944]
Iterations: 201
Function evaluations: 2171
Gradient evaluations: 201
Optimization Complete
-----------------------------------
Minimum found at (z1,z2,x) = (5.973519, 0.000000, 0.000000)
Coupling vars: 1.000000, 1.000000
Minimum objective: 1.36787944117
Elapsed time: 21578.5069466 milliseconds
What is it that I am doing wrong here? Also, are residuals defined in the component class reduced by Solver or Driver?
Upvotes: 1
Views: 169
Answers (1)
Kenneth Moore
Reputation: 2202
EDIT -- I initially implemented All-At-Once, so I've rewritten this. It should now be SAND.
So in my understanding of SAND, the optimizer minimizes the problem by varying the design variables simultaneously with the coupling variables to achieve feasibility and drive the residual constraint to zero.. This means the residual needs to be expressed explicitly so we don't need any implicit components or a solver -- the optimizer does it all. I modified your code as follows:
import time
import numpy as np
from openmdao.api import Component, Group, Problem, IndepVarComp, ExecComp, NLGaussSeidel, ScipyGMRES, \
ScipyOptimizer, pyOptSparseDriver
class SellarDis1(Component):
def __init__(self):
super(SellarDis1, self).__init__()
self.add_param('z', val=np.zeros(2))
self.add_param('x', val=0.0)
self.add_param('y2', val=1.0)
self.add_param('y1', val=1.0)
self.add_output('resid1', val=1.0)
def solve_nonlinear(self, params, unknowns, resids):
z1 = params['z'][0]
z2 = params['z'][1]
x1 = params['x']
y2 = params['y2']
y1 = params['y1']
unknowns['resid1'] = z1**2 + z2 + x1 - 0.2*y2 - y1
def linearize(self, params, unknowns, resids):
J = {}
J['resid1','y1'] = -1.0
J['resid1','y2'] = -0.2
J['resid1','z'] = np.array([[2*params['z'][0], 1.0]])
J['resid1','x'] = 1.0
return J
class SellarDis2(Component):
def __init__(self):
super(SellarDis2, self).__init__()
self.add_param('z', val=np.zeros(2))
self.add_param('y1', val=1.0)
self.add_param('y2', val=1.0)
self.add_output('resid2', val=1.0)
def solve_nonlinear(self, params, unknowns, resids):
z1 = params['z'][0]
z2 = params['z'][1]
y1 = params['y1']
y1 = abs(y1)
y2 = params['y2']
unknowns['resid2'] = y1**.5 + z1 + z2 - y2
def linearize(self, params, unknowns, resids):
J = {}
J['resid2', 'y2'] = -1.0
J['resid2', 'y1'] = 0.5*params['y1']**-0.5
J['resid2', 'z'] = np.array([[1.0, 1.0]])
return J
class SellarSAND(Group):
def __init__(self):
super(SellarSAND, self).__init__()
self.add('px', IndepVarComp('x', 1.0), promotes=['*'])
self.add('pz', IndepVarComp('z', np.array([5.0, 2.0])), promotes=['*'])
self.add('py1', IndepVarComp('y1', 1.0), promotes=['*'])
self.add('py2', IndepVarComp('y2', 1.0), promotes=['*'])
self.add('d1', SellarDis1(), promotes=['*'])
self.add('d2', SellarDis2(), promotes=['*'])
self.add('obj_cmp', ExecComp('obj = x**2 + z[1] + y1 + exp(-y2)',
z=np.array([0.0, 0.0]), x=0.0, y1=0.0, y2=0.0),
promotes=['*'])
self.add('con_cmp1', ExecComp('con1 = 3.16 - y1'), promotes=['*'])
self.add('con_cmp2', ExecComp('con2 = y2 - 24.0'), promotes=['*'])
top = Problem()
top.root = SellarSAND()
top.driver = ScipyOptimizer()
top.driver.options['optimizer'] = 'SLSQP'
top.driver.options['tol'] = 1.0e-12
top.driver.add_desvar('z', lower=np.array([-10.0, 0.0]),upper=np.array([10.0, 10.0]))
top.driver.add_desvar('x', lower=0.0, upper=10.0)
top.driver.add_desvar('y1', lower=-10.0, upper=10.0)
top.driver.add_desvar('y2', lower=-10.0, upper=10.0)
top.driver.add_objective('obj')
top.driver.add_constraint('con1', upper=0.0)
top.driver.add_constraint('con2', upper=0.0)
top.driver.add_constraint('resid1', equals=0.0)
top.driver.add_constraint('resid2', equals=0.0)
top.setup()
tt = time.time()
top.run()
print("\n")
print( "Minimum found at (z1,z2,x) = (%f, %f, %f)" % (top['z'][0], \
top['z'][1], \
top['x']))
print("Coupling vars: %f, %f" % (top['d1.y1'], top['d1.y2']))
print("Minimum objective: ", top['obj'])
print("Elapsed time: ", 1000*(time.time()-tt), "milliseconds")
You were on the right track; my main change was making it all explicit and declaring the y1 and y2 as des vars.
Doing this, I get
Minimum found at (z1,z2,x) = (1.977639, -0.000000, -0.000000)
Coupling vars: 3.160000, 3.755278
('Minimum objective: ', 3.1833939516406136)
('Elapsed time: ', 22.55988121032715, 'milliseconds')
Upvotes: 2
Related Questions
- ModuleNotFoundError using CaseReader in OpenMDAO
- Using APOPT (MINLP) solver in Openmdao
- Using pyOptSparseDriver with OpenMDAO and WISDEM
- Problem using ExternalCodeComp in Openmdao
- OpenMDAO efficiency with using multiple comp
- Openmdao 1.x: Efficient way to implement Expected Improvement
- Scaled paraboloid and derivatives checking
- Paraboloid optimization requiring scaling