этот вопрос соответствует этому один, но это не то же самое. Цель по-прежнему предназначена для студентов! Все еще играя с проблемой Селлара, я сравнил две разные проблемы:
Я ожидал получить аналогичные результаты с точки зрения обращений к дисциплинам (так как, по моим ожиданиям, если аналитические производные не даны и используется Ньютон, где-то нужно использовать FD для питания решателя Ньютона, и в этом случае, форсирование FD, когда даны производные, должно привести к аналогичному решению). Но задача 1 имеет следующее решение: Количество вызовов дисциплины 2: 9 а задача 2 имеет следующее решение: Количество вызовов дисциплин: 13
Поэтому обе проблемы не эквивалентны с точки зрения OpenMDAO. Это должно исходить из способа решения группы со связью с Ньютоном, когда аналитические производные не предоставляются, но я хотел бы понять, как это работает.
я не совсем понимаю, что вы делаете без производных? Определены ли аналитические производные или производные НЕ определены (вообще ничего не объявлено?)
Это было немного головокружительно. Ниже приведена автономная версия Sellar, которая работает на OpenMDAO V2.5, несмотря на использование Ньютон Солвер, в то время как НЕТ предоставляет какие-либо производные. Казалось бы, это вообще не должно работать, но работает (хотя и требует больше итераций, чем при объявлении производных с помощью 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)]
. Это означало, что вычисляемое обновление на каждой итерации было просто равно отрицательному значению остатка. Математически это фактически то же самое, что и решение с использованием решателя НелинейныйБлокЯкоби.
Это несколько неинтуитивное поведение произошло из-за того, что ExplicitComponent внутренне обрабатывает неявное преобразование и связанную с ним производную. Однако, если бы вместо этого вы определили компоненты Sellar как подклассы ImplicitComponent, тогда не объявление производных не сработало бы. Также обратите внимание, что вы не смогли бы оптимизировать эту модель без производных FD-d. В данном случае MDA заработало лишь из-за особенностей реализации ExplicitComponent.
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'])
Хорошо, большое спасибо за объяснение! Кристально чистый сейчас
Цель состояла в том, чтобы показать студентам, почему кодирование производных дисциплин (если возможно) имеет большое значение для эффективной оптимизации многодисциплинарных рабочих процессов. Вот почему мы тестируем все возможности (с FD/с производными) с GS/Newton, чтобы показать, что количество обращений к дисциплинам может кардинально измениться. Поэтому мы тестируем даже странные конфигурации, которые не должны встречаться при оптимизации в реальной жизни!
Если я добавлю self.declare_partials('', '', method='fd') для каждой дисциплины в задаче 1 (MDA Селлара без производной информации о дисциплинах с решателем Newton как NonlinearSolver ), я получаю вызовы 13 дисциплин, что теперь логично: обе задачи имеют одинаковый ответ