Openmdao V1.7鞍MDF

Question

我foound出与鞍问题的MDA一些奇怪的OpenMDAO的文档页面上（http://openmdao.readthedocs.io/en/1.7.3/usr-guide/tutorials/sellar.html）

如果我提取代码，只运行MDA（在学科添加计数器） ，我发现调用的次数是不同学科之间的差异（d1学科的d2的两倍），这是预期的。有人有答案吗？

这是效果

耦合瓦尔：25.588303，12.058488 学科1号和2次调用（10,5）

这里是代码

# For printing, use this import if you are running Python 2.x from __future__ import print_function 


import numpy as np 

from openmdao.api import Component from openmdao.api import ExecComp, IndepVarComp, Group, NLGaussSeidel, \ 
         ScipyGMRES 

class SellarDis1(Component): 
    """Component containing Discipline 1.""" 

    def __init__(self): 
     super(SellarDis1, self).__init__() 

     # Global Design Variable 
     self.add_param('z', val=np.zeros(2)) 

     # Local Design Variable 
     self.add_param('x', val=0.) 

     # Coupling parameter 
     self.add_param('y2', val=1.0) 

     # Coupling output 
     self.add_output('y1', val=1.0) 
     self.execution_count = 0 

    def solve_nonlinear(self, params, unknowns, resids): 
     """Evaluates the equation 
     y1 = z1**2 + z2 + x1 - 0.2*y2""" 

     z1 = params['z'][0] 
     z2 = params['z'][1] 
     x1 = params['x'] 
     y2 = params['y2'] 

     unknowns['y1'] = z1**2 + z2 + x1 - 0.2*y2 
     self.execution_count += 1 
    def linearize(self, params, unknowns, resids): 
     """ Jacobian for Sellar discipline 1.""" 
     J = {} 

     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): 
    """Component containing Discipline 2.""" 

    def __init__(self): 
     super(SellarDis2, self).__init__() 

     # Global Design Variable 
     self.add_param('z', val=np.zeros(2)) 

     # Coupling parameter 
     self.add_param('y1', val=1.0) 

     # Coupling output 
     self.add_output('y2', val=1.0) 
     self.execution_count = 0 
    def solve_nonlinear(self, params, unknowns, resids): 
     """Evaluates the equation 
     y2 = y1**(.5) + z1 + z2""" 

     z1 = params['z'][0] 
     z2 = params['z'][1] 
     y1 = params['y1'] 

     # Note: this may cause some issues. However, y1 is constrained to be 
     # above 3.16, so lets just let it converge, and the optimizer will 
     # throw it out 
     y1 = abs(y1) 

     unknowns['y2'] = y1**.5 + z1 + z2 
     self.execution_count += 1 
    def linearize(self, params, unknowns, resids): 
     """ Jacobian for Sellar discipline 2.""" 
     J = {} 

     J['y2', 'y1'] = .5*params['y1']**-.5 

     #Extra set of brackets below ensure we have a 2D array instead of a 1D array 
     # for the Jacobian; Note that Jacobian is 2D (num outputs x num inputs). 
     J['y2', 'z'] = np.array([[1.0, 1.0]]) 

     return J 



class SellarDerivatives(Group): 
    """ Group containing the Sellar MDA. This version uses the disciplines 
    with derivatives.""" 

    def __init__(self): 
     super(SellarDerivatives, self).__init__() 

     self.add('px', IndepVarComp('x', 1.0), promotes=['x']) 
     self.add('pz', IndepVarComp('z', np.array([5.0, 2.0])), promotes=['z']) 

     self.add('d1', SellarDis1(), promotes=['z', 'x', 'y1', 'y2']) 
     self.add('d2', SellarDis2(), promotes=['z', 'y1', 'y2']) 

     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=['obj', 'z', 'x', 'y1', 'y2']) 

     self.add('con_cmp1', ExecComp('con1 = 3.16 - y1'), promotes=['y1', 'con1']) 
     self.add('con_cmp2', ExecComp('con2 = y2 - 24.0'), promotes=['con2', 'y2']) 

     self.nl_solver = NLGaussSeidel() 
     self.nl_solver.options['atol'] = 1.0e-12 

     self.ln_solver = ScipyGMRES() 
     from openmdao.api import Problem, ScipyOptimizer 

top = Problem() top.root = SellarDerivatives() 

#top.driver = ScipyOptimizer() 
#top.driver.options['optimizer'] = 'SLSQP' 
#top.driver.options['tol'] = 1.0e-8 
# 
#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() 

# Setting initial values for design variables top['x'] = 1.0 top['z'] = np.array([5.0, 2.0]) 

top.run() 

print("\n") 

print("Coupling vars: %f, %f" % (top['y1'], top['y2'])) 


count1 = top.root.d1.execution_count 
count2 = top.root.d2.execution_count 
print("Number of discipline 1 and 2 calls (%i,%i)"% (count1,count2)) 

Answer 1

这是一个好的观察。每当你有一个循环，“头”组件再次运行。其原因如下：

如果您有包含隐性状态的分量的模型，一个执行是这样的：

1. 呼叫solve_nonlinear执行部件
2. 呼叫apply_nonlinear计算残差。

在这个模型中，我们没有任何隐式状态的组件，但是我们通过创建一个循环间接地创建了一个需要。我们的执行如下所示：

1. 致电solve_nonlinear执行所有组件。
2. 通过调用apply_nonlinear（将未知数缓存，请求solve_nolinear，并将差异保存在未知数中）仅放在“头部”组件上，以生成可以收敛的残差。

在这里，头部组件仅仅是第一个被执行的组件，然而它决定了循环运行的顺序。您可以通过构建一个更多的循环来验证只有一个头部组件获得额外的运行比2个组件。