有限差分OpenMDAO隐式状态

Question

我想用一个隐式变量和一个非线性求解器来迭代它的值，在一个系统中使用有限差分进行导数计算。我知道牛顿解需要解析导数，但是非线性高斯-西德尔似乎也没有正确的迭代。为了确保我已经适当地编写了代码，我在包中提供的intersect_parabola_line.py示例上测试了这种方法。

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上使用Python2.7.13和OpenMDAO 1.7.3。

Answer 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'