2017-04-25 1 views
1

私はコードを抽出し、唯一の(分野のカウンタを追加)MDAを実行した場合、私はOpenMDAOのドキュメントページ(http://openmdao.readthedocs.io/en/1.7.3/usr-guide/tutorials/sellar.htmlOpenmdao V1.7トルコ鞍MDF

にトルコ鞍問題のMDAとの奇妙な何かをfoound私は、呼び出しの数は、期待されていない分野間で異なることを観察します(d1分野のd2の2倍)。誰かが答えを持っていますか? 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)) 

答えて

1

これは、次のとおりです。ここで

は結果

カップリングはvarsのです良い観察。サイクルがあるたびに、「ヘッド」コンポーネントが2回実行されます。その理由は以下の通りである:

あなたが暗黙の状態が含まれているコンポーネントでモデルを持っている場合は、単一の実行は次のようになります。

  1. 残差を計算するためのコンポーネント
  2. コールapply_nonlinearを実行するためのコールsolve_nonlinear

このモデルでは暗黙の状態のコンポーネントはありませんが、サイクルを持つことによって間接的に必要なものが作成されました。実行は次のようになります。

  1. すべてのコンポーネントを実行するには、solve_nonlinearを呼び出します。
  2. (未知数をキャッシュし、solve_nolinearを呼び出し、未知数の差を保存する)を「ヘッド」コンポーネントに呼び出して、収束できる残差を生成します。ここで

、ヘッド部品は、それがサイクルを実行するために、どのような順序を決定しかしに基づいて実行されるだけで最初のコンポーネントである。あなたは、単一のヘッド部品が複数とのサイクルを構築することにより、余分な実行を取得していることを確認することができます2つのコンポーネントより。

+0

ありがとうございます!私は今行動をよく理解する。 – ThierryONERA

+0

私は今行動をよく理解しています。しかし、stateconnectionコンポーネントと一緒にNewtonアプローチでチュートリアルを使用すると、 'apply_nonlinear'関数を持つコンポーネントを追加したので、2つの分野は 'solve_linear'の部分に対してのみ呼び出されることになります。しかし、ここでも、規律1は2回と呼ばれます(シーケンスはdiscipline 1、discipline2、discipline 1、state connectionです)。 openmdaoにはまだ "ヘッド"コンポーネントが必要ですか? – ThierryONERA

+0

状態接続コンポーネントでは、まだモデルにサイクルがあると思います。状態接続コンポーネント(独自のapply_linearを持つ)が順序の "先頭"または最初のコンポーネントになるように、実行順序を手動で設定することによって、規律コンポーネントのいずれかが二重実行されないようにすることができます。 –

関連する問題