Ниже приведен пример кода, в котором BSplineComp сочетается либо с ExplicitComp, либо с ExternalCodeComp. Оба они выполняют одинаковые вычисления, и оба градиента компонентов рассчитываются с использованием конечной разности.
Если я запускаю версию Bspline+ExplicitComp, результат достигается за 2,3 итерации. Если я запускаю версию Bspline+ExternalCodeComp, мне приходится долго ждать. В этом случае он пытается найти градиент выхода по отношению к каждому входу. Так, например, есть 9 контрольных точек, которые интерполируются до 70 точек в компоненте bspline. Затем внешний компонент должен быть оценен столько, сколько интерполированных точек (70 раз)
Таким образом, в случае, когда bspline сочетается с дорогостоящим внешним кодом, для конечной разницы требуется столько же, сколько точек, в которые она интерполируется, что становится узким местом вычислений.
На основании этого ввода у меня есть два вопроса
1- Если компонент внешнего кода основан на явном компоненте, в чем основное отличие, вызывающее разницу в поведении? (учитывая, что оба имеют вход shape=70)
2- В ранее упомянутом сценарии, где bspline сочетается с дорогостоящим внешним кодом, был бы более эффективный способ их объединения, кроме того, как это показано здесь.
ГЛАВНЫЙ КОД: переменная «внешний» — это флаг для переключения внешнего/явного кода. установите значение true/false для запуска двух описанных выше случаев.
from openmdao.components.bsplines_comp import BsplinesComp
from openmdao.api import IndepVarComp, Problem, ExplicitComponent,ExecComp,ExternalCodeComp
from openmdao.api import ScipyOptimizeDriver, SqliteRecorder, CaseReader
import matplotlib.pyplot as plt
import numpy as np
external=True # change this to true for the case with external code comp. or false for the case with explicit comp.
rr=np.arange(0,70,1)
"Explicit component for the area under the line calculation"
class AreaComp(ExplicitComponent):
def initialize(self):
self.options.declare('lenrr', int)
self.options.declare('rr', types=np.ndarray)
def setup(self):
self.add_input('h', shape=lenrr)
self.add_output('area')
self.declare_partials(of='area', wrt='h', method='fd')
def compute(self, inputs, outputs):
rr = self.options['rr']
outputs['area'] = np.trapz(rr,inputs['h'])
class ExternalAreaComp(ExternalCodeComp):
def setup(self):
self.add_input('h', shape=70)
self.add_output('area')
self.input_file = 'paraboloid_input.dat'
self.output_file = 'paraboloid_output.dat'
# providing these is optional; the component will verify that any input
# files exist before execution and that the output files exist after.
self.options['external_input_files'] = [self.input_file]
self.options['external_output_files'] = [self.output_file]
self.options['command'] = [
'python', 'extcode_paraboloid.py', self.input_file, self.output_file
]
# this external code does not provide derivatives, use finite difference
self.declare_partials(of='*', wrt='*', method='fd')
def compute(self, inputs, outputs):
h = inputs['h']
# generate the input file for the paraboloid external code
np.savetxt(self.input_file,h)
# the parent compute function actually runs the external code
super(ExternalAreaComp, self).compute(inputs, outputs)
# parse the output file from the external code and set the value of f_xy
f_xy=np.load('a.npy')
outputs['area'] = f_xy
prob = Problem()
model = prob.model
n_cp = 9
lenrr = len(rr)
"Initialize the design variables"
x = np.random.rand(n_cp)
model.add_subsystem('px', IndepVarComp('x', val=x))
model.add_subsystem('interp', BsplinesComp(num_control_points=n_cp,
num_points=lenrr,
in_name='h_cp',
out_name='h'))
if external:
comp=ExternalAreaComp()
model.add_subsystem('AreaComp', comp)
else:
comp = AreaComp(lenrr=lenrr, rr=rr)
model.add_subsystem('AreaComp', comp)
case_recorder_filename2 = 'cases4.sql'
recorder2 = SqliteRecorder(case_recorder_filename2)
comp.add_recorder(recorder2)
comp.recording_options['record_outputs']=True
comp.recording_options['record_inputs']=True
model.connect('px.x', 'interp.h_cp')
model.connect('interp.h', 'AreaComp.h')
model.add_constraint('interp.h', lower=0.9, upper=1, indices=[0])
prob.driver = ScipyOptimizeDriver()
prob.driver.options['optimizer'] = 'SLSQP'
prob.driver.options['disp'] = True
#prob.driver.options['optimizer'] = 'COBYLA'
#prob.driver.options['disp'] = True
prob.driver.options['tol'] = 1e-9
model.add_design_var('px.x', lower=1,upper=10)
model.add_objective('AreaComp.area',scaler=1)
prob.setup(check=True)
#prob.run_model()
prob.run_driver()
cr = CaseReader(case_recorder_filename2)
case_keys = cr.system_cases.list_cases()
cou=-1
for case_key in case_keys:
cou=cou+1
case = cr.system_cases.get_case(case_key)
plt.plot(rr,case.inputs['h'],'-*')
Внешний код extcode_paraboloid.py ниже
import numpy as np
if __name__ == '__main__':
import sys
input_filename = sys.argv[1]
output_filename = sys.argv[2]
h=np.loadtxt(input_filename)
rr=np.arange(0,70,1)
rk= np.trapz(rr,h)
np.save('a',np.array(rk))