To give you a bit of a background, I am trying to maximize the total coverage by a satellite constellation. So what I basically do is estimate the initial state vector for every satellite, propagate them over a period of time and co-relate the area covered with my area of interest and estimate the coverage. The initial state vector is an array of 6 elements (3 position + 3 velocity). So if I have 2 satellites then the array will have 2 rows and 6 columns. The below simple test script mimics the same data flow and formulation of my actual code to calculate the initial position.
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'])
This code works as intended. The shape of d(phase angle) and rv0 depend on Nd. Now to optimize this I will have to consider Nd (number of satellites) as a design variable (i.e declare it as ivc and have an other component or execcomp to calculate d). This however gives me an error as the shape of parameters cannot change after setup() . Is there a workaround for this ? I could calculate the position (and with it the coverage) for every single satellite by doing the model and problem setup inside a for loop and combine them in the end. But then again how can I access Nd (for the for loop) and make coverage the objective.
Edit: Now that I think about this, I do know the maximum value Nd can be and I could use this define the shape of the parameters. For example if max is 10 the I would define rv0 as self.add_output('rv0', val=np.zeros([10,6])). So if Nd is 2 in a iteration then outputs['rv0'][:2,:] = rv0.real. The rest will be zero. I dont think this will have any impact on the overall coverage. But from Openmdao point will this cause any issue, that I dont see ? Or Is there a better way of doing this ?
There are two basic options:
I would personally chose option 2 for a couple of reasons. I would think the relaxation based approach would be more effective for larger numbers of N (MINLP are very hard to solve well). Also, its the most flexible approach, since you don't need to re-allocate anything. You can just adjust the values of the active array accordingly and turn things on an off at will.
Here is some updated code where you can do either method:
import numpy as np
class InitialState(om.ExplicitComponent):
def initialize(self):
self.options.declare('max_Nd', default=10, types=int,
desc='maximum number of allowed satellites')
def setup(self):
Nd = self.options['max_Nd']
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(Nd))
self.add_input('active', val=np.zeros(Nd)) # range from 0 to 1 to activate or de-activate
self.add_output('rv0', val=np.zeros([Nd,6]))
#self.declare_partials(of='*', wrt='*', method='fd')
self.declare_partials('*', '*', method='cs')
def compute(self, inputs, outputs):
Nd = self.options['max_Nd']
a = inputs['a']
b = inputs['b']
c = inputs['c']
d = inputs['d']
active = inputs['active']
for i in range(Nd):
outputs['rv0'][i,:3] = [a**2,b*c*d[i],0]
outputs['rv0'][i,3:] = [d[i]**2,a*b,0]
outputs['rv0'][i,:] *= active[i]
model = om.Group()
Nd = 4
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))
ivc.add_output('active', val=np.zeros(Nd))
model.add_subsystem('Pos', InitialState(max_Nd = Nd))
model.connect('a', 'Pos.a')
model.connect('b', 'Pos.b')
model.connect('c', 'Pos.c')
model.connect('d', 'Pos.d')
model.connect('active', 'Pos.active')
prob= om.Problem(model)
prob.setup()
prob['d'] = np.random.randint(10, size=(Nd))
prob['active'] = np.round(np.random.random(Nd))
prob.run_model()
print(prob['Pos.rv0'])