2017-11-09 7 views
0

微分計算に有限差分を使用するシステムで、暗黙変数と非線形ソルバーを使用してその値を反復したいとします。私はニュートンソルバーが分析微分を必要としていることを理解していますが、非線形Gauss-Siedelは正しく反復していないようです。適切にコード化したことを確かめるために、このアプローチはパッケージに含まれているintersect_parabola_line.pyの例でテストしました。私が手有限差分のOpenMDAO暗黙状態

from __future__ import print_function 
from openmdao.api import Component, Group, Problem, Newton, ScipyGMRES,NLGaussSeidel 

class Line(Component): 
    """Evaluates y = -2x + 4.""" 
    def __init__(self): 
     super(Line, self).__init__() 
     self.add_param('x', 1.0) 
     self.add_output('y', 0.0) 
     # User can change these. 
     self.slope = -2.0 
     self.intercept = 4.0 

    def solve_nonlinear(self, params, unknowns, resids): 
     """ y = -2x + 4 """ 
     x = params['x'] 
     m = self.slope 
     b = self.intercept 
     unknowns['y'] = m*x + b 

class Parabola(Component): 
    """Evaluates y = 3x^2 - 5""" 
    def __init__(self): 
     super(Parabola, self).__init__() 
     self.add_param('x', 1.0) 
     self.add_output('y', 0.0) 
     # User can change these. 
     self.a = 3.0 
     self.b = 0.0 
     self.c = -5.0 

    def solve_nonlinear(self, params, unknowns, resids): 
     """ y = 3x^2 - 5 """ 
     x = params['x'] 
     a = self.a 
     b = self.b 
     c = self.c 
     unknowns['y'] = a*x**2 + b*x + c 

class Balance(Component): 
    """Evaluates the residual y1-y2""" 
    def __init__(self): 
     super(Balance, self).__init__() 
     self.add_param('y1', 0.0) 
     self.add_param('y2', 0.0) 
     self.add_state('x', 5.0) 

    def solve_nonlinear(self, params, unknowns, resids): 
     """This component does no calculation on its own. It mainly holds the 
     initial value of the state. An OpenMDAO solver outside of this 
     component varies it to drive the residual to zero.""" 
     pass 

    def apply_nonlinear(self, params, unknowns, resids): 
     """ Report the residual y1-y2 """ 
     y1 = params['y1'] 
     y2 = params['y2'] 
     resids['x'] = y1 - y2 

if __name__ == '__main__': 
    top = Problem() 
    root = top.root = Group() 
    root.add('line', Line()) 
    root.add('parabola', Parabola()) 
    root.add('bal', Balance()) 

    root.connect('line.y', 'bal.y1') 
    root.connect('parabola.y', 'bal.y2') 
    root.connect('bal.x', 'line.x') 
    root.connect('bal.x', 'parabola.x') 
    root.deriv_options['type'] = 'fd' 

    root.nl_solver = NLGaussSeidel() #Newton() 
    root.ln_solver = ScipyGMRES() 
    root.nl_solver.options['iprint'] = 2 

    top.setup() 

    # Positive solution 
    top['bal.x'] = 7.0 
    root.list_states() 
    top.run() 
    print('Positive Solution x=%f, line.y=%f, parabola.y=%f' % (top['bal.x'], top['line.y'], top['parabola.y'])) 

    # Negative solution 
    top['bal.x'] = -7.0 
    root.list_states() 
    top.run() 
    print('Negative Solution x=%f, line.y=%f, parabola.y=%f' % (top['bal.x'], top['line.y'], top['parabola.y'])) 

出力は次のようになります。

States in model: 

bal.x 
Value: 7.0 
Residual: 0.0 

[root] NL: NLN_GS 1 | 152 1 
[root] NL: NLN_GS 2 | 152 1 
[root] NL: NLN_GS 2 | Converged in 2 iterations 
Positive Solution x=7.000000, line.y=-10.000000, parabola.y=142.000000 

States in model: 

bal.x 
Value: -7.0 
Residual: -152.0 

[root] NL: NLN_GS 1 | 124 1 
[root] NL: NLN_GS 2 | 124 1 
[root] NL: NLN_GS 2 | Converged in 2 iterations 
Negative Solution x=-7.000000, line.y=18.000000, parabola.y=142.000000 

任意のヒントをいただければ幸いです。私はOSXでPython 2.7.13とOpenMDAO 1.7.3を使用しています。

答えて

1

したがって、非線形Gauss Seidelは、モデルをバランス成分と暗黙の状態で実際に収束させることはできません。あなたは本当にニュートンを使ってください。その作業をするには、ニュートンが収束するモデルも有限差分を使用することを確認する必要があります。ルートグループにfdを設定すると、そのシステムよりも上にある微分計算では、そのシステム全体で近似されたfdが表示されます。

は、この作業を行うには、これを試してみてください。

top = Problem() 
root = top.root = Group() 
comp1 = root.add('line', Line()) 
comp2 = root.add('parabola', Parabola()) 
comp3 = root.add('bal', Balance()) 

root.connect('line.y', 'bal.y1') 
root.connect('parabola.y', 'bal.y2') 
root.connect('bal.x', 'line.x') 
root.connect('bal.x', 'parabola.x') 
root.deriv_options['type'] = 'fd' 
comp1.deriv_options['type'] = 'fd' 
comp2.deriv_options['type'] = 'fd' 
comp3.deriv_options['type'] = 'fd' 
+0

私は 'root.deriv_options'は、すべてのコンポーネントに流れ落ちたと思っていたが、明らかに私が間違っていました。このヒントをありがとうございます。私の実装をもっときれいにするでしょう。乾杯。 – gbarter

関連する問題