Это было немного головокружение наверняка.Ниже приведена автономная версия sellar, которая работает на OpenMDAO V2.5, несмотря на использование NewtonSolver , тогда как NOT предоставляет любые производные.Это, по-видимому, не должно работать вообще, но это работает (хотя это и требует больше итераций, чем когда вы объявляли производные с FD) !!.
Так что здесь происходит немного неуловимо и является функциейо том, как ExplicitComponent на самом деле реализован в рамках OpenMDAO.Я отошлю вас к статье OpenMDAO для получения более подробной информации, но 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'])