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個組件。