Чтобы дать вам некоторую предысторию, я пытаюсь максимизировать общий охват спутниковой группировкой. Итак, что я в основном делаю, так это оцениваю начальный вектор состояния для каждого спутника, распространяю их в течение определенного периода времени и связываю область, покрытую с моей интересующей областью, и оцениваю покрытие. Вектор начального состояния представляет собой массив из 6 элементов (3 позиции + 3 скорости). Итак, если у меня есть 2 спутника, тогда в массиве будет 2 строки и 6 столбцов. Приведенный ниже простой тестовый сценарий имитирует тот же поток данных и формулировку моего фактического кода для вычисления начальной позиции.
import openmdao.api as om
import numpy as np
class InitialState(om.ExplicitComponent):
def __init__(self, Nd):
super(InitialState, self).__init__()
self.Nd = Nd
def setup(self):
self.add_input('a', val=120.456)
self.add_input('b', val=58)
self.add_input('c', val=46)
self.add_input('d', val=np.zeros(self.Nd))
self.add_output('rv0', val=np.zeros([self.Nd,6]))
#self.declare_partials(of='*', wrt='*', method='fd')
self.declare_partials('*', '*')
def compute(self, inputs, outputs):
d = inputs['d']
rv0 = np.zeros([len(d),6])
for i in range(len(d)):
r0, v0 = self.calc(inputs['a'], inputs['b'],
inputs['c'], d[i])
rv0[i,:3] = r0
rv0[i,3:] = v0
outputs['rv0'] = rv0.real
def calc(self, a, b, c, d):
r0 = np.array([a**2,b*c*d,0],dtype=complex)
v0 = np.array([d**2,a*b,0],dtype=complex)
return r0,v0
def compute_partials(self, inputs, J):
h = 1e-16
ih = complex(0, h)
rv_drn = np.zeros(4,complex)
Jacs = np.zeros((6, 4))
for j in range(len(inputs['d'])):
rv_drn[:] = [inputs['a'], inputs['b'], inputs['c'], inputs['d'][j]]
start = j*6
stop = start+6
for i in range(4):
rv_drn[i] += ih
r0, v0 = self.calc(rv_drn[0], rv_drn[1], rv_drn[2], rv_drn[3])
rv_drn[i] -= ih
Jacs[:3, i] = r0.imag/h
Jacs[3:, i] = v0.imag/h
J['rv0', 'a'][start:stop] = [[k] for k in Jacs[:, 0]]
J['rv0', 'b'][start:stop] = [[k] for k in Jacs[:, 1]]
J['rv0', 'c'][start:stop] = [[k] for k in Jacs[:, 2]]
J['rv0', 'd'][start:stop] = [[k] for k in Jacs[:, 3]]
model = om.Group()
Nd = 2
ivc = model.add_subsystem('ivc', om.IndepVarComp(), promotes_outputs=['*'])
ivc.add_output('a', val=120.456)
ivc.add_output('b', val=58)
ivc.add_output('c', val=46)
ivc.add_output('d', val=np.zeros(Nd))
model.add_subsystem('Pos', InitialState(Nd))
model.connect('a', 'Pos.a')
model.connect('b', 'Pos.b')
model.connect('c', 'Pos.c')
model.connect('d', 'Pos.d')
prob= om.Problem(model)
prob.setup()
prob['d'] = np.random.randint(10, size=(Nd))
prob.check_partials(compact_print=True)
prob.run_model()
print(prob['Pos.rv0'])
Этот код работает так, как задумано. Форма d (фазовый угол) и rv0 зависит от Nd. Теперь, чтобы оптимизировать это, мне нужно будет рассматривать Nd (количество спутников) как переменную проекта (т.е. объявить его как iv c и иметь другой компонент или execcomp для вычисления d). Однако это дает мне ошибку, поскольку форма параметров не может измениться после setup(). Есть ли обходной путь для этого? Я мог рассчитать положение (а вместе с ним и покрытие) для каждого отдельного спутника, выполнив настройку модели и задачи внутри для l oop и комбинируя их в конце. Но опять же, как я могу получить доступ к Nd (для for loop) и сделать охват целью.
Изменить: теперь, когда я думаю об этом, я знаю, что максимальное значение Nd может быть и Я мог бы использовать это для определения формы параметров. Например, если max равно 10, я бы определил rv0 как self.add_output('rv0', val=np.zeros([10,6])). Итак, если Nd равно 2 на итерации, то outputs['rv0'][:2,:] = rv0.real. Остальное будет ноль. Я не думаю, что это повлияет на общий охват. Но с точки зрения Openmdao это вызовет какие-либо проблемы, которых я не вижу? Или есть лучший способ сделать это?