(Openmdao 2.4.0) 组件的 'compute_partials' 功能似乎在为该组件强制 'declare_partials' 到 FD 时运行

Question

我想使用 Group 的 Newton 非线性求解器为 Sellar 求解 MDA。我已经定义了带有导数的学科（使用'compute_partials'），但我想检查在强制学科不使用其分析导数（使用问题中的'declare_partials'时）调用学科'compute'和'compute_partials'的次数定义 ）。问题是，即使我强制不使用它，似乎仍然调用了“compute_partials”函数。 这是一个例子（Sellar）

所以对于纪律 2，我添加了一个计数器，我有

from openmdao.test_suite.components.sellar import SellarDis1, SellarDis2 

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

    def _do_declares(self):
        # Analytic Derivs
        self.declare_partials(of='*', wrt='*')
        self.exec_count_d = 0

    def compute_partials(self, inputs, J):
        """
        Jacobian for Sellar discipline 2.
        """
        y1 = inputs['y1']
        if y1.real < 0.0:
            y1 *= -1

        J['y2', 'y1'] = .5*y1**-.5
        J['y2', 'z'] = np.array([[1.0, 1.0]])
        self.exec_count_d += 1

我创建了一个与 on OpendMDAO docs 类似的 MDA，但调用我创建的 SellarDis2withDerivatives 和 SellarDis1withDerivatives 并像这样更改 Newton_solver() 的非线性求解器

    cycle.add_subsystem('d1', SellarDis1withDerivatives(), promotes_inputs=['x', 'z', 'y2'], promotes_outputs=['y1'])
    cycle.add_subsystem('d2', SellarDis2withDerivatives(), promotes_inputs=['z', 'y1'], promotes_outputs=['y2'])

    # Nonlinear Block Gauss Seidel is a gradient free solver
    cycle.nonlinear_solver = NewtonSolver()
    cycle.linear_solver = DirectSolver()

然后我运行以下问题

 prob2 = Problem()

prob2.model = SellarMDA()

prob2.setup()

prob2.model.cycle.d1.declare_partials('*', '*', method='fd')
prob2.model.cycle.d2.declare_partials('*', '*', method='fd')

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

prob2.run_model()


count = prob2.model.cycle.d2.exec_count_d
print("Number of derivatives calls (%i)"% (count))

结果，我得到了

=====

周期

NL：牛顿在 3 次迭代中收敛 衍生品调用次数 (3)

因此，似乎仍然以某种方式调用了函数“compute_partials”（即使导数是使用 FD 计算的）。有人作为解释吗？

Answer 1

我认为这是一个错误（或者可能是指定导数的意外结果。）

这种行为是导数混合声明的副产品，我们允许用户将组件上的某些导数指定为“fd”，而将其他导数指定为解析。因此，我们总是能够在一个组件上同时执行 fd 和 compute_partials。

我们可以在 openmdao 中进行两个更改来解决这个问题：

1. 如果没有明确声明为解析的导数，请勿调用 compute_partials。

2. 过滤掉任何声明为 'fd' 的变量，这样如果用户尝试在 compute_partials 中设置它们，就会引发 keyerror（或者可能只是警告，并且不会覆盖导数值）

与此同时，唯一的解决方法是注释掉 compute_partials 方法，或者将组件包含在一个组中并对该组进行有限差分。

Answer 2

另一种解决方法是在你的类中有一个属性（这里称为_call_compute_partials），它跟踪是否声明了任何分析导数。并且compute_partials() 中的条件可以在调用方法的方法之外实现。

from openmdao.core.explicitcomponent import ExplicitComponent
from openmdao.core.indepvarcomp import IndepVarComp
from openmdao.core.problem import Problem
from openmdao.drivers.scipy_optimizer import ScipyOptimizeDriver


class ExplicitComponent2(ExplicitComponent):

    def __init__(self, **kwargs):
        super(ExplicitComponent2, self).__init__(**kwargs)
        self._call_compute_partials = False

    def declare_partials(self, of, wrt, dependent=True, rows=None, cols=None, val=None,
                         method='exact', step=None, form=None, step_calc=None):
        if method == 'exact':
            self._call_compute_partials = True
        super(ExplicitComponent2, self).declare_partials(of, wrt, dependent, rows, cols, val,
                         method, step, form, step_calc)


class Cylinder(ExplicitComponent2):
    """Main class"""

    def setup(self):
        self.add_input('radius', val=1.0)
        self.add_input('height', val=1.0)

        self.add_output('Area', val=1.0)
        self.add_output('Volume', val=1.0)

        # self.declare_partials('*', '*', method='fd')
        # self.declare_partials('*', '*')

        self.declare_partials('Volume', 'height', method='fd')
        self.declare_partials('Volume', 'radius', method='fd')
        self.declare_partials('Area', 'height', method='fd')
        self.declare_partials('Area', 'radius')
        # self.declare_partials('Area', 'radius', method='fd')

    def compute(self, inputs, outputs):
        radius = inputs['radius']
        height = inputs['height']

        area = height * radius * 2 * 3.14 + 3.14 * radius ** 2 * 2
        volume = 3.14 * radius ** 2 * height
        outputs['Area'] = area
        outputs['Volume'] = volume

    def compute_partials(self, inputs, partials):
        if self._call_compute_partials:
            print('Calculate partials...')


if __name__ == "__main__":

    prob = Problem()

    indeps = prob.model.add_subsystem('indeps', IndepVarComp(), promotes=['*'])
    indeps.add_output('radius', 2.)  # height
    indeps.add_output('height', 3.)  # radius
    main = prob.model.add_subsystem('cylinder', Cylinder(), promotes=['*'])

    # setup the optimization
    prob.driver = ScipyOptimizeDriver()

    prob.model.add_design_var('radius', lower=0.5, upper=5.)
    prob.model.add_design_var('height', lower=0.5, upper=5.)
    prob.model.add_objective('Area')
    prob.model.add_constraint('Volume', lower=10.)

    prob.setup()
    prob.run_driver()
    print(prob['Volume'])  # should be around 10