这个问题与这个一是一致的,但是不一样。目标仍然是为了学生的目的!我还在玩塞拉问题,我比较了两个不同的问题:
我希望在学科调用方面得到类似的结果(因为,在我的预期中,如果没有给出解析导数,并且在某个地方使用牛顿,那么就必须用FD来给牛顿求解器,在这种情况下,当给出导数时强迫FD会导致类似的解)。但问题1有以下解决方案: discipline2调用数:9;问题2有以下解决方案:学科调用数: 13
因此,从OpenMDAO的角度来看,这两个问题并不等价。当不提供解析导数时,它应该来自于与牛顿耦合的群的求解方法,但我想了解它是如何工作的。
发布于 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工作起来。
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'])
https://stackoverflow.com/questions/54686634
复制相似问题