（Openmdao 2.4.0）在衍生品学科上不提供衍生品/强制FD的区别

Question

这个问题符合这个one，但它不一样。目标仍然是为了学生的目的！仍在玩塞拉尔问题，我比较了两个不同的问题：

• 问题1：Sellar的MDA没有使用Newton求解器作为NonlinearSolver的Disciplines的衍生信息
• 问题2：Sellar的MDA，带有牛顿求解器的学科的衍生信息为NonlinearSolver，但在问题级别的每个学科上都有选项declare_partials（''，''，method ='fd'）
• 对于两者，linearsolver是相同的，并且都计算类似点的MDA

我希望在对学科的调用方面获得类似的结果（正如我所期望的那样，如果没有给出分析导数并且使用牛顿，在某个地方，必须使用FD来提供牛顿求解器，在这种情况下，在给出衍生物时强制FD应该导致类似的解决方案）。但问题1有以下解决方案：纪律2呼叫次数：9和问题2有以下解决方案：纪律调用次数：13

因此，从OpenMDAO的角度来看，这两个问题都不相同。当没有提供分析衍生物时，它应该来自解决与牛顿耦合的群体的方式，但我想了解它是如何工作的。

Answer 1

这肯定是一个令人头疼的问题。下面是一个自包含的sellar版本，适用于OpenMDAO V2.5，尽管使用NewtonSolver而不提供任何衍生产品。这似乎根本不应该工作，但确实如此（虽然它确实需要比你用FD声明衍生物时更多的迭代）！

所以这里发生的事情有点微妙，并且是在OpenMDAO的封面下实际实现ExplicitComponent的函数。我将向您介绍OpenMDAO paper以获取更多详细信息，但OpenMDAO实际上将所有内容转换为隐式表单。所以每个显式输出实际上都得到R(output_var) = compute(inputs)[output_var] - outputs[output_var]形式的残差。

那么当你运行牛顿时会发生什么，调用计算函数，然后形成残差驱动输出变量向量以匹配计算值。这是使用标准的Newton方法完成的：[dR/du] [delta-u] = -[R(u)]。

那么，如果你不提供任何衍生工具，它是如何工作的呢？好吧，请注意dR_i/du_i = -1（这是显式变量的残差的导数，相对于输出向量中的关联值）。 OpenMDAO ExplicitComponent类自动定义这个偏导数。关于输入的其他衍生物随后由ExplicitComponent的子类提供。因此，当你没有定义任何衍生物时，你仍然得到了dR_i/du_i = -1。

然后牛顿法退化为-[I] [delta-u] = -[R(u)]。这意味着每次迭代的计算更新仅等于残差的负值。在数学上，这实际上与使用NonlinearBlockJacobi求解器求解相同。

这种有点不直观的行为发生是因为ExplicitComponent在内部处理隐式转换和相关的派生本身。但是，如果您已将Sellar组件定义为ImplicitComponent的子类，则不会声明衍生系统不起作用。另外，请注意，如果没有FD-d衍生物，您将无法使用此模型进行优化。这只是ExplicitComponent实现的一个怪癖，它使MDA在这种情况下工作。

import numpy as np

from openmdao.api import ExplicitComponent, Group, Problem, NewtonSolver, \
                         DirectSolver, IndepVarComp, ExecComp

class SellarDis1(ExplicitComponent):
    """
    Component containing Discipline 1 -- no derivatives version.
    """

    def setup(self):

        # Global Design Variable
        self.add_input('z', val=np.zeros(2))

        # Local Design Variable
        self.add_input('x', val=0.)

        # Coupling parameter
        self.add_input('y2', val=1.0)

        # Coupling output
        self.add_output('y1', val=1.0)

        # Finite difference all partials.
        # self.declare_partials('*', '*', method='fd')

    def compute(self, inputs, outputs):
        """
        Evaluates the equation
        y1 = z1**2 + z2 + x1 - 0.2*y2
        """
        z1 = inputs['z'][0]
        z2 = inputs['z'][1]
        x1 = inputs['x']
        y2 = inputs['y2']

        outputs['y1'] = z1**2 + z2 + x1 - 0.2*y2
        print('compute y1', outputs['y1'])

class SellarDis2(ExplicitComponent):
    """
    Component containing Discipline 2 -- no derivatives version.
    """

    def setup(self):
        # Global Design Variable
        self.add_input('z', val=np.zeros(2))

        # Coupling parameter
        self.add_input('y1', val=1.0)

        # Coupling output
        self.add_output('y2', val=1.0)

        # Finite difference all partials.
        # self.declare_partials('*', '*', method='fd')

    def compute(self, inputs, outputs):
        """
        Evaluates the equation
        y2 = y1**(.5) + z1 + z2
        """

        z1 = inputs['z'][0]
        z2 = inputs['z'][1]
        y1 = inputs['y1']
        print('y1', 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
        if y1.real < 0.0:
            y1 *= -1

        outputs['y2'] = y1**.5 + z1 + z2

class SellarMDA(Group):
    """
    Group containing the Sellar MDA.
    """

    def setup(self):
        indeps = self.add_subsystem('indeps', IndepVarComp(), promotes=['*'])
        indeps.add_output('x', 1.0)
        indeps.add_output('z', np.array([5.0, 2.0]))

        cycle = self.add_subsystem('cycle', Group(), promotes=['*'])
        cycle.add_subsystem('d1', SellarDis1(), promotes_inputs=['x', 'z', 'y2'], promotes_outputs=['y1'])
        cycle.add_subsystem('d2', SellarDis2(), promotes_inputs=['z', 'y1'], promotes_outputs=['y2'])

        # Nonlinear Block Gauss Seidel is a gradient free solver
        newton = cycle.nonlinear_solver = NewtonSolver()
        newton.options['iprint'] = 2
        newton.options['maxiter'] = 20
        newton.options['solve_subsystems'] = False
        cycle.linear_solver = DirectSolver()

        self.add_subsystem('obj_cmp', ExecComp('obj = x**2 + z[1] + y1 + exp(-y2)',
                           z=np.array([0.0, 0.0]), x=0.0),
                           promotes=['x', 'z', 'y1', 'y2', 'obj'])

        self.add_subsystem('con_cmp1', ExecComp('con1 = 3.16 - y1'), promotes=['con1', 'y1'])
        self.add_subsystem('con_cmp2', ExecComp('con2 = y2 - 24.0'), promotes=['con2', 'y2'])


prob = Problem()

prob.model = SellarMDA()

prob.setup()

prob['x'] = 2.
prob['z'] = [-1., -1.]

prob.run_model()


print(prob['y1'])
print(prob['y2'])