(Openmdao 2.4.0) разница между предоставлением производных / принудительным FD для дисциплин с производными

этот вопрос соответствует этому один, но это не то же самое. Цель по-прежнему предназначена для студентов! Все еще играя с проблемой Селлара, я сравнил две разные проблемы:

  • Задача 1: MDA Селлара без производной информации о дисциплинах с решателем Ньютона как NonlinearSolver
  • Проблема 2: MDA Селлара с производной информацией о дисциплинах с решателем Ньютона как NonlinearSolver, но с опциями declare_partials('', '', method='fd') для каждой дисциплины на уровне задачи
  • для обоих линейный решатель один и тот же, и оба вычисляют MDA для аналогичной точки

Я ожидал получить аналогичные результаты с точки зрения обращений к дисциплинам (так как, по моим ожиданиям, если аналитические производные не даны и используется Ньютон, где-то нужно использовать FD для питания решателя Ньютона, и в этом случае, форсирование FD, когда даны производные, должно привести к аналогичному решению). Но задача 1 имеет следующее решение: Количество вызовов дисциплины 2: 9 а задача 2 имеет следующее решение: Количество вызовов дисциплин: 13

Поэтому обе проблемы не эквивалентны с точки зрения OpenMDAO. Это должно исходить из способа решения группы со связью с Ньютоном, когда аналитические производные не предоставляются, но я хотел бы понять, как это работает.

Если я добавлю self.declare_partials('', '', method='fd') для каждой дисциплины в задаче 1 (MDA Селлара без производной информации о дисциплинах с решателем Newton как NonlinearSolver ), я получаю вызовы 13 дисциплин, что теперь логично: обе задачи имеют одинаковый ответ

ThierryONERA 14.02.2019 11:38

я не совсем понимаю, что вы делаете без производных? Определены ли аналитические производные или производные НЕ определены (вообще ничего не объявлено?)

Justin Gray 14.02.2019 13:03
Стоит ли изучать PHP в 2023-2024 годах?
Стоит ли изучать PHP в 2023-2024 годах?
Привет всем, сегодня я хочу высказать свои соображения по поводу вопроса, который я уже много раз получал в своем сообществе: "Стоит ли изучать PHP в...
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
В JavaScript одним из самых запутанных понятий является поведение ключевого слова "this" в стрелочной и обычной функциях.
Приемы CSS-макетирования - floats и Flexbox
Приемы CSS-макетирования - floats и Flexbox
Здравствуйте, друзья-студенты! Готовы совершенствовать свои навыки веб-дизайна? Сегодня в нашем путешествии мы рассмотрим приемы CSS-верстки - в...
Тестирование функциональных ngrx-эффектов в Angular 16 с помощью Jest
В системе управления состояниями ngrx, совместимой с Angular 16, появились функциональные эффекты. Это здорово и делает код определенно легче для...
Концепция локализации и ее применение в приложениях React ⚡️
Концепция локализации и ее применение в приложениях React ⚡️
Локализация - это процесс адаптации приложения к различным языкам и культурным требованиям. Это позволяет пользователям получить опыт, соответствующий...
Пользовательский скаляр GraphQL
Пользовательский скаляр GraphQL
Листовые узлы системы типов GraphQL называются скалярами. Достигнув скалярного типа, невозможно спуститься дальше по иерархии типов. Скалярный тип...
2
2
82
1

Ответы 1

Это было немного головокружительно. Ниже приведена автономная версия 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'])

Хорошо, большое спасибо за объяснение! Кристально чистый сейчас

ThierryONERA 15.02.2019 14:32

Цель состояла в том, чтобы показать студентам, почему кодирование производных дисциплин (если возможно) имеет большое значение для эффективной оптимизации многодисциплинарных рабочих процессов. Вот почему мы тестируем все возможности (с FD/с производными) с GS/Newton, чтобы показать, что количество обращений к дисциплинам может кардинально измениться. Поэтому мы тестируем даже странные конфигурации, которые не должны встречаться при оптимизации в реальной жизни!

ThierryONERA 15.02.2019 14:49

Другие вопросы по теме