具有有限差分的 OpenMDAO 隐式状态

Question

我想使用隐式变量和非线性求解器在使用有限差分进行导数计算的系统中迭代其值。我知道牛顿求解器需要解析导数，但非线性 Gauss-Siedel 似乎也没有正确迭代。为了确保我已经正确编码，我在包中提供的 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 上使用 Python 2.7.13 和 OpenMDAO 1.7.3。

Answer 1

因此，非线性高斯赛德尔并不能真正收敛具有平衡分量和隐式状态的模型。你真的应该使用牛顿。为了实现这一点，您还需要确保牛顿收敛的模型也使用有限差分。当您在根组中设置 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'