首页
学习
活动
专区
圈层
工具
发布
首页
学习
活动
专区
圈层
工具
MCP广场
社区首页 >问答首页 >(Openmao2.4.0)在带导数的学科上不提供导数/强迫FD之间的区别

(Openmao2.4.0)在带导数的学科上不提供导数/强迫FD之间的区别
EN

Stack Overflow用户
提问于 2019-02-14 09:03:39
回答 1查看 89关注 0票数 2

这个问题与这个是一致的,但是不一样。目标仍然是为了学生的目的!我还在玩塞拉问题,我比较了两个不同的问题:

  • 问题1:以牛顿求解器为NonlinearSolver的不含导数信息的鞍区模型
  • 问题2:以牛顿求解器为NonlinearSolver,但在问题级别上每个学科的选项declare_partials('','',method='fd')上具有导数信息的问题2。
  • 对于这两种情况,线性求解器都是相同的,并且都计算出类似点的MDA值。

我希望在学科调用方面得到类似的结果(因为,在我的预期中,如果没有给出解析导数,并且在某个地方使用牛顿,那么就必须用FD来给牛顿求解器,在这种情况下,当给出导数时强迫FD会导致类似的解)。但问题1有以下解决方案: discipline2调用数:9;问题2有以下解决方案:学科调用数: 13

因此,从OpenMDAO的角度来看,这两个问题并不等价。当不提供解析导数时,它应该来自于与牛顿耦合的群的求解方法,但我想了解它是如何工作的。

EN

回答 1

Stack Overflow用户

发布于 2019-02-14 15:23:42

这肯定有点伤脑筋。下面是一个在OpenMDAO V2.5上工作的自成体系的塞勒版本,尽管它使用的是NewtonSolver,而不是提供了任何衍生品。这似乎根本不应该起作用,但它确实起作用(尽管它比使用FD声明衍生产品时需要更多的迭代时间)!

所以,这里发生的事情有点微妙,这是一个函数,它是如何在ExplicitComponent的掩护下在OpenMDAO中实际实现的。我将向OpenMDAO纸提供更多细节,但OpenMDAO实际上将所有内容转换为隐藏的隐式形式。因此,每个显式输出实际上都得到了表单R(output_var) = compute(inputs)[output_var] - outputs[output_var]的残差。

所以当你运行牛顿时,会调用计算函数,然后形成一个残差,驱动输出变量向量与计算值相匹配。这是用标准牛顿法:[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内部处理隐式转换和相关的导数本身。然而,如果你把塞勒组件定义为ImplicitComponent的子类,那么不声明衍生工具就行不通了。另外,请注意,如果没有FD-d导数,您就无法用这个模型进行优化。在这种情况下,这只是ExplicitComponent实现的一个怪癖,使得MDA工作起来。

代码语言:javascript
运行
复制
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'])
票数 1
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/54686634

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档