Я хочу решить MDA для Селлара, используя нелинейный решатель Newton для группы. Я определил Дисциплины с производными (используя 'compute_partials'), но я хочу проверить количество вызовов дисциплины 'compute' и 'compute_partials' при принуждении или отказе дисциплин от использования их аналитических производных (используя 'declare_partials' в задаче определение ). Проблема в том, что кажется, что функция 'compute_partials' все еще вызывается, хотя я заставляю ее не использовать. Вот пример (Селлар)
Итак, для Дисциплины 2 я добавляю счетчик, и у меня есть
from openmdao.test_suite.components.sellar import SellarDis1, SellarDis2
class SellarDis2withDerivatives(SellarDis2):
"""
Component containing Discipline 2 -- derivatives version.
"""
def _do_declares(self):
# Analytic Derivs
self.declare_partials(of='*', wrt='*')
self.exec_count_d = 0
def compute_partials(self, inputs, J):
"""
Jacobian for Sellar discipline 2.
"""
y1 = inputs['y1']
if y1.real < 0.0:
y1 *= -1
J['y2', 'y1'] = .5*y1**-.5
J['y2', 'z'] = np.array([[1.0, 1.0]])
self.exec_count_d += 1
Я создаю аналогичный MDA как в документах OpendMDAO, но вызываю SellarDis2withDerivatives, который я создал, и SellarDis1withDerivatives и меняю linear_solver для Newton_solver() следующим образом
cycle.add_subsystem('d1', SellarDis1withDerivatives(), promotes_inputs=['x', 'z', 'y2'], promotes_outputs=['y1'])
cycle.add_subsystem('d2', SellarDis2withDerivatives(), promotes_inputs=['z', 'y1'], promotes_outputs=['y2'])
# Nonlinear Block Gauss Seidel is a gradient free solver
cycle.nonlinear_solver = NewtonSolver()
cycle.linear_solver = DirectSolver()
Затем я запускаю следующую проблему
prob2 = Problem()
prob2.model = SellarMDA()
prob2.setup()
prob2.model.cycle.d1.declare_partials('*', '*', method='fd')
prob2.model.cycle.d2.declare_partials('*', '*', method='fd')
prob2['x'] = 2.
prob2['z'] = [-1., -1.]
prob2.run_model()
count = prob2.model.cycle.d2.exec_count_d
print("Number of derivatives calls (%i)"% (count))
И в результате получаю
=====
NL: Ньютон сошелся за 3 итерации Количество вызовов деривативов (3)
Поэтому кажется, что функция 'compute_partials' все равно как-то вызывается (даже если производные вычисляются с помощью FD). Кто-нибудь в качестве объяснения?
Я считаю, что это ошибка (или, возможно, непреднамеренное следствие того, как определяются производные).
Такое поведение является побочным продуктом смешанного объявления производных, когда мы позволяем пользователю указывать некоторые производные для компонента как «fd», а другие производные — как аналитические. Таким образом, мы всегда можем выполнить как fd, так и compute_partials
на компоненте.
Чтобы исправить это, мы можем внести два изменения в openmdao:
Не вызывайте compute_partials
, если никакие производные не были явно объявлены аналитическими.
Отфильтруйте любые переменные, объявленные как «fd», чтобы при попытке пользователя установить их в compute_partials
возникала ошибка ключа (или, может быть, просто предупреждение, а производное значение не перезаписывалось)
В то же время единственным обходным решением будет закомментировать метод compute_partials
или, альтернативно, заключить компонент в группу, а конечную разность - в группу.
Спасибо. Я собирался спросить, для какого варианта использования потребуется возможность конечной разности компонента, для которого определены аналитические производные, вне контекста проверки партиалов. Так что это скорее показательный случай.
Другим обходным решением является наличие атрибута (здесь он называется _call_compute_partials
) в вашем классе, который отслеживает, если там объявлены какие-либо аналитические производные. И условное выражение в compute_partials()
может быть реализовано вне метода, где метод вызывается.
from openmdao.core.explicitcomponent import ExplicitComponent
from openmdao.core.indepvarcomp import IndepVarComp
from openmdao.core.problem import Problem
from openmdao.drivers.scipy_optimizer import ScipyOptimizeDriver
class ExplicitComponent2(ExplicitComponent):
def __init__(self, **kwargs):
super(ExplicitComponent2, self).__init__(**kwargs)
self._call_compute_partials = False
def declare_partials(self, of, wrt, dependent=True, rows=None, cols=None, val=None,
method='exact', step=None, form=None, step_calc=None):
if method == 'exact':
self._call_compute_partials = True
super(ExplicitComponent2, self).declare_partials(of, wrt, dependent, rows, cols, val,
method, step, form, step_calc)
class Cylinder(ExplicitComponent2):
"""Main class"""
def setup(self):
self.add_input('radius', val=1.0)
self.add_input('height', val=1.0)
self.add_output('Area', val=1.0)
self.add_output('Volume', val=1.0)
# self.declare_partials('*', '*', method='fd')
# self.declare_partials('*', '*')
self.declare_partials('Volume', 'height', method='fd')
self.declare_partials('Volume', 'radius', method='fd')
self.declare_partials('Area', 'height', method='fd')
self.declare_partials('Area', 'radius')
# self.declare_partials('Area', 'radius', method='fd')
def compute(self, inputs, outputs):
radius = inputs['radius']
height = inputs['height']
area = height * radius * 2 * 3.14 + 3.14 * radius ** 2 * 2
volume = 3.14 * radius ** 2 * height
outputs['Area'] = area
outputs['Volume'] = volume
def compute_partials(self, inputs, partials):
if self._call_compute_partials:
print('Calculate partials...')
if __name__ == "__main__":
prob = Problem()
indeps = prob.model.add_subsystem('indeps', IndepVarComp(), promotes=['*'])
indeps.add_output('radius', 2.) # height
indeps.add_output('height', 3.) # radius
main = prob.model.add_subsystem('cylinder', Cylinder(), promotes=['*'])
# setup the optimization
prob.driver = ScipyOptimizeDriver()
prob.model.add_design_var('radius', lower=0.5, upper=5.)
prob.model.add_design_var('height', lower=0.5, upper=5.)
prob.model.add_objective('Area')
prob.model.add_constraint('Volume', lower=10.)
prob.setup()
prob.run_driver()
print(prob['Volume']) # should be around 10
Хорошо, я лучше понимаю, как это работает. Вопрос был немного сложным, потому что цель состояла в том, чтобы объяснить студентам влияние Ньютона / Гаусса Зейделя с производными / без них, но начиная с того же Компонента и MDA, поэтому просто изменяя выбор в определении проблемы.