如何使用 OpenMDAO 实现 SAND 架构
How to implement SAND architecture using OpenMDAO
我正在尝试在 OpenMDAO 上使用 Sellar Problem 实施同步分析和设计 (SAND) 架构。我想到了以下方法 -
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")
但是在 运行 系统上我得到的结果不正确 -
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
我这里做错了什么？
此外，组件 class 中定义的残差是否被 Solver 或 Driver 减少？
编辑——我最初实现了一次全部，所以我重写了它。现在应该是 SAND。
所以在我对 SAND 的理解中，优化器通过同时改变设计变量和耦合变量来最小化问题，以实现可行性并将残差约束驱动为零。这意味着残差需要明确表达，所以我们不需要任何隐式组件或求解器——优化器可以完成所有工作。我修改了你的代码如下：
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")
你走在正确的轨道上；我的主要更改是使其全部明确并将 y1 和 y2 声明为 des vars。
这样做，我得到
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')
我正在尝试在 OpenMDAO 上使用 Sellar Problem 实施同步分析和设计 (SAND) 架构。我想到了以下方法 -
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")
但是在 运行 系统上我得到的结果不正确 -
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
我这里做错了什么？ 此外，组件 class 中定义的残差是否被 Solver 或 Driver 减少？
编辑——我最初实现了一次全部，所以我重写了它。现在应该是 SAND。
所以在我对 SAND 的理解中，优化器通过同时改变设计变量和耦合变量来最小化问题，以实现可行性并将残差约束驱动为零。这意味着残差需要明确表达，所以我们不需要任何隐式组件或求解器——优化器可以完成所有工作。我修改了你的代码如下：
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")
你走在正确的轨道上；我的主要更改是使其全部明确并将 y1 和 y2 声明为 des vars。
这样做，我得到
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')