如何使用OpenMDAO实现SAND体系结构

我正试图通过Sellar问题在OpenMDAO上实现同步分析和设计（SAND）体系结构。我想下面的方式这样做的 -如何使用OpenMDAO实现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

它是什么，我做错了什么？另外，由Solver或Driver减少的组件类中定义的残差？

来源

2016-02-05 Divya Manglam

编辑 - 我最初实施一次完成，所以我重写了这一点。它现在应该是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')

来源

2016-02-05 14:56:36

注，编辑我的反应，使代码是沙而不是AAO。 –

谢谢。顺便说一下，对于这个问题，AAO和SAND有什么不同的实施方式。 AAO应具有实际和目标耦合变量以及残差约束之间的兼容性约束。但是对于优化器来说，这不会是多余的吗？ –

对于AAO，我制作了单独的组件来容纳残差，保留在d1和2之间的数据传输中，并且学科d1和d2的代码与MDF示例中的代码相同。所以，差异是微妙的，但多存储一点。总体而言，我认为约束的数量是相同的。 –

如何使用OpenMDAO实现SAND体系结构

回答

相关问题