输入参数数组的可变大小

Question

为了给您一些背景知识，我正在尝试最大化卫星星座的总覆盖范围。所以我基本上做的是估计每颗卫星的初始状态向量，在一段时间内传播它们并将覆盖的区域与我感兴趣的区域相关联并估计覆盖范围。初始状态向量是一个包含 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（卫星数量）作为设计变量（即将其声明为 ivc 并使用其他组件或 execcomp 来计算 d）。然而，这给了我一个错误，因为参数的形状在 setup() 之后不能改变。有解决方法吗？我可以通过在 for 循环中进行模型和问题设置并最终将它们组合起来，来计算每颗卫星的位置（以及覆盖范围）。但是话又说回来，我怎样才能访问Nd（对于for loop）并覆盖objective。

编辑：现在我想到了这个，我知道最大值 Nd 可以，我可以用它来定义参数的形状。例如，如果最大值为 10，我会将 rv0 定义为 self.add_output('rv0', val=np.zeros([10,6]))。因此，如果 Nd 在一次迭代中为 2，则 outputs['rv0'][:2,:] = rv0.real。其余的将为零。我认为这不会对整体覆盖范围产生任何影响。但是从 Openmdao 的角度来看，这会导致任何我看不到的问题吗？或者有更好的方法吗？

Answer 1

有两个基本选项：

将卫星数量作为您组的参数，并分配正确大小的数组。然后运行对每个 N 值进行单独优化，然后选择您最喜欢的一个。
为您要考虑的尽可能多的卫星分配阵列，并添加一个额外的消隐阵列作为输入以打开和关闭某些卫星。然后您可以使用混合整数方法或尝试一些基于松弛的方法并将消隐数组视为实数。

出于几个原因，我个人会选择选项 2。我认为基于松弛的方法对于更大数量的 N 会更有效（MINLP 很难很好地解决）。此外，它是最灵活的方法，因为您不需要重新分配任何东西。您可以相应地调整 active 数组的值并随意打开或关闭。

这是一些更新的代码，您可以在其中执行任一方法：

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'])

输入参数数组的可变大小

Variable size for input parameter array

openmdao