我的设计变量对OpenMDAO中的约束和目标函数没有影响

问题描述 投票：0回答：1

所以我正在研究 OpenMDAO 的优化问题。优化问题的主题是根据您选择的某些任务要求找到最佳体积，然后使用两种计算权重和计算增量的方法使体积收敛。所以最后我的问题是在完成代码后，我运行它，然后收到一条错误消息：“DerivativesWarning：约束或目标 [('obj.obj', inds=[0]), ('c1. c1', inds=[0]), ('c2.c2', inds=[0])] 不会受到影响 由问题的设计变量决定。”和“DerivativesWarning：设计变量 [('V', inds=[0]), ('FR', inds=[0])] 对约束或目标没有影响。”

所以在我的代码中有 3 个子类计算 2 个重量和 1 个关于空气动力学值。然后在我的优化类中，我放置了 3 个子系统，连接我的输入和输出，然后创建我的目标函数和约束，最后我通过实现优化来解决问题，然后我尝试运行它并收到错误消息 所以这是代码（它有很多公式，但我认为这对我的主题来说并不重要）：

import openmdao.api as om
from math import *

AR_tails = 1.0
pho_cruise_alt = 0.00211 #in slug/ft3 4000ft
pho_max_alt = 0.001876 #in slug/ft3 8000ft
pho_sea_level = 0.002377 #in slug/ft3
Vs = 84.5 #en f/s
mu = 3.66*10**(-7) #in poise
tc = 0.15 #no unit
nb_lobes = 3 #no unit a voir si on veut mettre ca dasn l'optimisation
range = 1000 #in nautical miles 
NE = 4 #number of engine, no unit, a voir par la suite avec la reliaility
BR = 0.7 #no unit
BSFC = 0.48 #Brake Specific fuel consumption, no unit
Payload = 40000 #in lb
FS = 4 #Security factor, no unit
Manufacturing_factor = 1.2 #no unit
nb_septum = 2 #no unit,a voir si je modifie ca pour permettre de le modif 
Faf = 1.26 #account for the external attach fitting, no unit
Seaming_factor = 1.06 #no unit
nb_ball = nb_lobes*2 #6 #no unit
Fpsq = 1.0 #account for the fabrication of the tail in lightweight structural material, no unit
surface_CS = 0.2 #percentage of the CS on the tail surface, no unit
Fact = 0.79 #no unit
Finstall = 1.15 #account for the installation of actuators, no unit
lec = 150 #account for the length of control cable for the engines, in ft
Nb_prop = 4 #no unit
Nbl = 3 #no unit
Kp = 31.92 #no unit
nt = 2 #number of fuel tank, no unit
Integral_tank = 0 #percentage of fuel tank that are integral, no unit
Npad = 3 #number of pads for ACLS, no unit
Vsr = 4 #landing sink rate, in f/s
Wcomputer = 250 #account for the weight of the computer, it's estimated, in lb
Wavionics = 250 #account for the weight of instruments, in lb
Wfs = 470 #account for the weight of fuel system, in lb
Kelect = 33.73 #account for a coef that depend on the type of range flight, no unit
np = 0.65 #propeller efficiency, no unit


class Geometry(om.ExplicitComponent):
    def setup(self):
        self.add_input('V',val=1000000.0)
        self.add_input('FR',val=3.0)
        self.add_output('V23',val=1.0)
        self.add_output('de',val=3.0)
        self.add_output('lb',val=1.0)
        self.add_output('dc',val=1.0)
        self.add_output('w',val=1.0)
        self.add_output('ht',val=1.0)
        self.add_output('AR',val=1.0)
        self.add_output('Swet_body_lobed',val=1.0)
        self.add_output('SHT',val=1.0)
        self.add_output('SVT',val=1.0)
        self.add_output('Wg',val=1.0)
        self.add_output('Wh0',val=1.0)
        self.add_output('Wh1',val=1.0)
        self.add_output('total_fuel',val=1.0)
        self.add_output('Woe',val=1.0)
        self.add_output('Cd0',val=1.0)
        self.add_output('K',val=1.0)

    def setup_partials(self):
        # Finite difference all partials.
        self.declare_partials('*', '*', method='fd')

    def compute(self,inputs,outputs):
        V = inputs['V']
        FR = inputs['FR']
        V23 = outputs['V23']
        de = outputs['de']
        lb = outputs['lb']
        dc = outputs['dc']
        w = outputs['w']
        ht = outputs['ht']
        AR = outputs['AR']
        Swet_body_lobed = outputs['Swet_body_lobed']
        SHT = outputs['SHT']
        SVT = outputs['SVT']
        Wg= outputs['Wg']
        Wh0=outputs['Wh0']
        Wh1=outputs['Wh1']
        total_fuel=outputs['total_fuel']
        Woe = outputs['Woe']
        Cd0 = outputs['Cd0']
        K = outputs['K']

        V23 = (float(V))**(2/3)
        de = ((6*V)/(FR*pi))**(1/3) #equivalent diameter of body of revolution
        lb = FR*de
        dc = (de/(-0.0178*nb_lobes**2+0.361*nb_lobes+0.575))
        w = (1+nb_lobes)*(dc/2)
        ht = dc
        AR = ((4*w**2)/(pi*lb*w))
        p = 1.6075
        Swet_body_ellipsoid = pi*((lb**p*w**p+lb**p*ht**p+w**p*ht**p)/3)**(1/p)
        if nb_lobes == 1:
            p_ellipsoid = 2*pi*(dc/2)
            p_lobes = 2*pi*(dc/2)
        if nb_lobes == 2:
            p_ellipsoid = 2.525*pi*(dc/2)
            p_lobes = (8*pi*(dc/2))/3
        if nb_lobes == 3:
            p_ellipsoid = 3.084*pi*(dc/2)
            p_lobes = (10*pi*(dc/2))/3
        if nb_lobes == 4:
            p_ellipsoid = 3.663*pi*(dc/2)
            p_lobes = (12*pi*(dc/2))/3
        if nb_lobes == 5:
            p_ellipsoid = 4.254*pi*(dc/2)
            p_lobes = (14*pi*(dc/2))/3
        Swet_body_lobed = (p_lobes/p_ellipsoid)*Swet_body_ellipsoid
        CHT = -0.0051*(1000000/V)+0.0717
        SHT = CHT*(V23*lb)/(0.38*lb)
        CVT = -0.0049*(1000000/V)+0.0641
        SVT = CVT*(V23*lb)/(0.38*lb)
        Cd0tab = []
        q = (1/2)*pho_cruise_alt*Vs**2
        Re = ((lb*Vs*pho_cruise_alt)/mu)
        Cf = (0.455/(log10(Re)**2.58))
        FF3D_body = 1+(1.5/FR**1.5)+(7/FR**3)
        Cd0_body = ((FF3D_body*Cf*Swet_body_lobed)/V23)
        Cd0tab.append(Cd0_body)
        FF_tail = 1+1.2*tc+100*tc**4
        ctail = ((AR_tails*SHT*0.5)**0.5+(AR_tails*SVT*0.5)**0.5)/2
        Re_tail = (ctail*Vs*pho_cruise_alt)/mu
        cf_tail = 0.455/(log10(Re_tail)**2.58)
        Swet_tail = 2.2*(SHT+SVT)
        Cd0_tail = (FF_tail*cf_tail*Swet_tail)/V23
        Cd0tab.append(Cd0_tail)
        #Add every Cd0 for each part of the airship, check what we choose in the end 
        Cd0cab = (0.108*Cd0_body*V23+7.7)/V23
        Cd0tab.append(Cd0cab)
        Cd0eng = NE*4.25/V23
        Cd0tab.append(Cd0eng)
        Cd0engcooling = NE*(2*10**(-6)*V+4.1)/V23
        Cd0tab.append(Cd0engcooling)
        Cd0engmount = (0.044*Cd0_body*V23+0.92)/V23
        Cd0tab.append(Cd0engmount)
        Cd0cables = (9.7*10**(-6)*V+10.22)/V23
        Cd0tab.append(Cd0cables)
        Cd0ACLS = 0.0002
        Cd0tab.append(Cd0ACLS)
        Cd0_int = (4.78*10**(-6)*V)/V23
        Cd0tab.append(Cd0_int)
        Cd0 = sum(Cd0tab)
        K = -0.0145*(1/AR)**4 + 0.182*(1/AR)**3 - 0.514*(1/AR)**2 + 0.838*(1/AR) - 0.053
        if nb_lobes == 1:
            Nl = 2
        if nb_lobes == 2:
            Nl = 2.25
        if nb_lobes == 3:
            Nl = 2.4
        if nb_lobes == 4:
            Nl = 2.5
        if nb_lobes == 5:
            Nl = 2.54
        K = K/Nl
        L = 0.0646*V*(pho_max_alt/pho_sea_level)
        Wlanding = L/BR
        Wh1 = Wlanding - L
        A = (326*np)/(BSFC*sqrt(K*Cd0))
        B = q*V23*sqrt(Cd0/K)
        Wh0 = B*tan((range/A)+atan(Wh1/B))
        fuel_burned = Wh0-Wh1
        fuelres = fuel_burned*0.05
        total_fuel = fuel_burned+fuelres
        Wzf = (L/BR)-fuelres
        Woe = Wzf - Payload
        Wg = Wlanding + fuel_burned
       

class aerodynamic(om.ExplicitComponent):
    def setup(self):
        self.add_input('Wh0',val=1.0)
        self.add_input('V23',val=1.0)
        self.add_input('Cd0',val=1.0)
        self.add_input('K',val=1.0)
        self.add_output('qmax',val=1.0)
        self.add_output('Maxpower_perengine',val=1.0)
        self.add_output('Dp',val=1.0)
        self.add_output('np_check',val=1.0)

    def setup_partials(self):
        self.declare_partials('*','*',method='fd')

    def compute(self,inputs,outputs):   
        Wh0= inputs['Wh0']
        V23= inputs['V23']
        Cd0= inputs['Cd0']
        K= inputs['K']
        qmax= outputs['qmax']
        Maxpower_perengine= outputs['Maxpower_perengine']
        Dp= outputs['Dp']
        np_check= outputs['np_check']


        L_aero = Wh0
        qmax = pho_sea_level*(((Vs*1.1)**2)/2)
        Clmax_power = (L_aero/qmax/V23)
        Drag = (Cd0+K*Clmax_power**2)*qmax*V23
        Maxpower_perengine = (Vs*1.1)*(Drag/np/NE/550)
        Cs = ((pho_sea_level*Vs**5)/Maxpower_perengine/550/10**2)**(1/5) #10 correspond a une rotation max pour le propeller speed and pho is in slug/ft3
        J = 0.156*Cs**2+0.241*Cs+0.138
        Dp = ((Vs)/10)/J
        np_check = 0.139*Cs**3-0.749*Cs**2+1.37*Cs+0.0115

class second_weight_estimation(om.ExplicitComponent):
    def setup(self):
        self.add_input('qmax',val=1.0)
        self.add_input('ht',val=1.0)
        self.add_input('V',val=1000000.0)
        self.add_input('SHT',val=1.0)
        self.add_input('SVT',val=1.0)
        self.add_input('Swet_body_lobed',val=1.0)
        self.add_input('dc',val=1.0)
        self.add_input('lb',val=1.0)
        self.add_input('Maxpower_perengine',val=1.0)
        self.add_input('Dp',val=1.0)
        self.add_input('total_fuel',val=1.0)
        self.add_input('Woe',val=1.0)
        self.add_input('Wh0',val=1.0)
        self.add_input('Wh1',val=1.0)
        self.add_output('Wg2',val=1.0)

    def setup_partials(self):
        self.declare_partials('*','*',method='fd')

    def compute(self,inputs,outputs):
        qmax = inputs['qmax']
        ht = inputs['ht']
        V = inputs['V']
        SHT = inputs['SHT']
        SVT = inputs['SVT']
        Swet_body_lobed = inputs['Swet_body_lobed']
        dc = inputs['dc']
        lb = inputs['lb']
        Maxpower_perengine = inputs['Maxpower_perengine']
        Dp = inputs['Dp']
        total_fuel = inputs['total_fuel']
        Woe = inputs['Woe']
        Wh0 = inputs['Wh0']
        Wh1 = inputs['Wh1']
        Wg2 = outputs['Wg2']



        Weight = []
        Pint = 1.2*qmax+0.0635*ht
        hull_fabric_load = FS*12*(Pint/144)*(dc/2)
        hull_fabric_density = 0.0085*hull_fabric_load+1.365
        Whull = hull_fabric_density*Manufacturing_factor*Faf*Swet_body_lobed/(16*9) #to convert each unit in yd2 and then go into lb
        Septum_fabric_load = 1.5*hull_fabric_load
        Septum_fabric_density = 0.0085*Septum_fabric_load+1.365
        side_body_area = 0.75 #percentage of envelope side body area 
        Wseptum = nb_septum*Seaming_factor*Septum_fabric_density*side_body_area*pi*ht*(lb/4)/(16*9) #we have a formula for the body side area after side body area
        Wenv = Wseptum+Whull
        Weight.append(Wenv)
        Vball = V*((1/(pho_max_alt/pho_sea_level))-1)
        Sball = nb_ball*pi*(nb_lobes*((Vball/pi)/6))**(2/3)
        Wball = 0.035*Sball
        Weight.append(Wball)
        Wssf = Fpsq*(SHT+SVT)*Faf*(1-surface_CS)
        Wcs = Fpsq*(SHT+SVT)*surface_CS
        Wtail = Wcs+Wssf
        Weight.append(Wtail)
        Wact = Fact*Finstall*surface_CS*(SHT+SVT)
        Weng = NE*Maxpower_perengine**(0.7956)*4.848
        Weight.append(Weng)
        Wengmt = 0.64*Weng
        Weight.append(Wengmt)
        Wec = 60.27*(lec*(NE/100))**(0.724)
        Weight.append(Wec)
        Wstart = 50.38*((Weng/1000))**(0.459)
        Weight.append(Wstart)
        Wprop = Nb_prop*Kp*Nbl**(0.391)*(Dp*(Maxpower_perengine/1000))**(0.782)
        Weight.append(Wprop)
        Wft = 2.49*(total_fuel/6)**(0.6)*nt**(0.2)*NE**(0.13)*(1/(1+Integral_tank))**(0.3) #total fuel is divided by 6 which the weight in lb per gallon for aviation gas
        Weight.append(Wft)
        Wpress_syst = 0.02*Woe
        Weight.append(Wpress_syst)
        Ppad = Pint
        ACLS_area_landing = 3*(0.23*Wh1*(Vsr/(Npad*Ppad))) #we multiply it by 3 because 3 pads but the configuration can change some stuff and then modify the calculations
        ACLS_area_takeoff = Wh0/Ppad
        Wacls = 1.6*max(ACLS_area_landing,ACLS_area_takeoff)
        Weight.append(Wacls)
        Wvms = Wcomputer+Wavionics+Wact
        Weight.append(Wvms)
        Welect = Kelect * (Wfs+Wcomputer+Wavionics)**(0.51)  #maybe don't have this equation 
        Weight.append(Welect)
        Wmsys = 0.05*Woe
        Weight.append(Wmsys)
        Wuf = 0.01*total_fuel
        Weight.append(Wuf)
        Wmargin = 0.06*Woe
        Weight.append(Wmargin)
        Woe2 = sum(Weight)
        Wg2 = Woe2 + total_fuel + Payload

class Optimization(om.Group):
    def setup(self):
        cycle = self.add_subsystem('cycle',om.Group(),promotes_inputs=['V','FR'])
        cycle.add_subsystem('d1', Geometry(),promotes_inputs=['V','FR'])
        cycle.add_subsystem('d4', aerodynamic())
        cycle.add_subsystem('d5', second_weight_estimation(),promotes_inputs=['V'])

        cycle.connect('d1.V23','d4.V23')
        cycle.connect('d1.lb','d5.lb')
        cycle.connect('d1.dc','d5.dc')
        cycle.connect('d1.ht','d5.ht')
        cycle.connect('d1.Swet_body_lobed','d5.Swet_body_lobed')
        cycle.connect('d1.SHT','d5.SHT')
        cycle.connect('d1.SVT','d5.SVT')
        cycle.connect('d1.Wh0','d4.Wh0')
        cycle.connect('d1.Wh0','d5.Wh0')
        cycle.connect('d1.Wh1','d5.Wh1')
        cycle.connect('d1.total_fuel','d5.total_fuel')
        cycle.connect('d1.Woe','d5.Woe')
        cycle.connect('d1.Cd0','d4.Cd0')
        cycle.connect('d1.K','d4.K')
        cycle.connect('d4.qmax','d5.qmax')
        cycle.connect('d4.Maxpower_perengine','d5.Maxpower_perengine')
        cycle.connect('d4.Dp','d5.Dp')


        cycle.nonlinear_solver = om.NonlinearBlockGS()

        cycle.set_input_defaults('V',val=2)

        self.add_subsystem('obj', om.ExecComp('obj = Wg - Wg2'), promotes= ['obj','Wg','Wg2'])
        self.add_subsystem('c1', om.ExecComp('c1 = L/Wg'), promotes=['c1','L','Wg'])
        self.add_subsystem('c2', om.ExecComp('c2 = np_check'), promotes=['c2','np_check'])


prob = om.Problem()
prob.model = Optimization()

prob.driver = om.ScipyOptimizeDriver()
prob.driver.options['optimizer'] = 'SLSQP'
# prob.driver.options['maxiter'] = 100
prob.driver.options['tol'] = 1e-8

prob.model.add_design_var('V', lower=0)
prob.model.add_design_var('FR', lower=0)
prob.model.add_objective('obj')
prob.model.add_constraint('c1', lower= 0, upper=0.6)
prob.model.add_constraint('c2', lower=0.95*np, upper=1.05*np)

# Ask OpenMDAO to finite-difference across the model to compute the gradients for the optimizer
prob.model.approx_totals()

prob.setup()
prob.set_val('V',2000000)
prob.set_val('FR',3)
prob.set_solver_print(level=0)

prob.run_driver()

print('minimum found at')
print(prob.get_val('V'))
print(prob.get_val('FR'))

print('minumum objective')
print(prob.get_val('obj'))

总错误消息：

C:\Users\umroot\anaconda3\envs\fastuav\lib\site-packages\openmdao\core\total_jac.py:1626: DerivativesWarning:Constraints or objectives [('obj.obj', inds=[0]), ('c1.c1', inds=[0]), ('c2.c2', inds=[0])] cannot be impacted 
by the design variables of the problem.
C:\Users\umroot\anaconda3\envs\fastuav\lib\site-packages\openmdao\core\total_jac.py:1655: DerivativesWarning:Design variables [('V', inds=[0]), ('FR', inds=[0])] have no impact on the constraints or objective.
Positive directional derivative for linesearch    (Exit mode 8)
            Current function value: 0.0
            Iterations: 5
            Function evaluations: 1
            Gradient evaluations: 1
Optimization FAILED.
Positive directional derivative for linesearch
-----------------------------------
minimum found at
[2000000.]
[3.]
minumum objective
[0.]

预先感谢您的查看并希望我们能够共同找到解决方案

optimization openmdao
1个回答
0
投票

发布的代码中存在一些问题。

当你写作时

np_check= outputs['np_check']

...

np_check = 5

Python 不会像您想象的那样将 5 写入 

outputs['np_check']
，而是将其重新定义为新整数。因此，outputs['np_check'] 保持其默认值（它在计算方法中永远不会被覆盖）。

要解决此问题，您可以直接分配给 

outputs['np_check']
，或者使用以下命令强制 numpy 将 5 插入 np_check 的数组值中：

np_check[...] = 5

这明确告诉 numpy 您正在将 5 推入 np_check 数组的所有元素中（它仅包含实例中的单个元素）。

在此之后，您仍然会收到一条警告，指出约束和目标不会受到设计变量的影响。在这种情况下，最好的办法是打开 OpenMDAO 在执行脚本的目录的 

reports/n2.html
中自动为您创建的 n 方图。

在此图中，变量列在左侧。任何橙色变量都不会连接到模型中其他位置的输出。对于 

L
和 
FR
这是有道理的，因为它们最终是优化器提供的设计变量。

但该图表明 

L
、
Wg
、
Wg2
和 
np_check
未连接到输出。 OpenMDAO 假设这些输入仅保留分配给它们的任何默认值（通常为 1.0），但优化器不会更改它们。

这些值似乎是 

cycle
组内的输出，但它们需要提升，首先从其计算组件中取出，然后从 
cycle
组本身中取出。

            cycle = self.add_subsystem('cycle',om.Group(),
                                       promotes_inputs=['V','FR'],
                                       promotes_outputs=['Wg', 'Wg2', 'L', 'np_check'])
            cycle.add_subsystem('d1', Geometry(),
                                promotes_inputs=['V','FR'],
                                promotes_outputs=['Wg', 'L'])
            cycle.add_subsystem('d4', aerodynamic(),
                                promotes_outputs=['np_check'])
            cycle.add_subsystem('d5', second_weight_estimation(),
                                promotes_inputs=['V'],
                                promotes_outputs=['Wg2'])

仅供参考，我通过单击 N2 图左侧菜单栏中的“circled i”工具发现了这一点，您可以看到模型中变量的“升级”名称。共享相同升级名称的变量在 OpenMDAO 中隐式连接。

我做的另一件事是从 

cycle
组中删除非线性求解器分配。该组没有反馈，因此首选默认非线性求解器。

最后，其中一个约束和目标的值似乎约为 1E6，因此我通过将这些响应的 

ref
值设置为 1.0E6 来缩放它们。这将导致优化器将 
1.0E6
的值视为 
1.0
。 
ref0
和 
add_objective
有一个对应的 
add_constraint
参数，它提供优化器视为零的值。

通过这些更改，优化器强制设计进入 

log10(Re)
无效的区域。发生这种情况时，您可以尝试引发 OpenMDAO AnalysisError，这将迫使某些优化器走回有效区域，但在这种情况下，它似乎坚持行驶到该区域。

    import openmdao.api as om
    from math import *

    AR_tails = 1.0
    pho_cruise_alt = 0.00211 #in slug/ft3 4000ft
    pho_max_alt = 0.001876 #in slug/ft3 8000ft
    pho_sea_level = 0.002377 #in slug/ft3
    Vs = 84.5 #en f/s
    mu = 3.66*10**(-7) #in poise
    tc = 0.15 #no unit
    nb_lobes = 3 #no unit a voir si on veut mettre ca dasn l'optimisation
    range = 1000 #in nautical miles
    NE = 4 #number of engine, no unit, a voir par la suite avec la reliaility
    BR = 0.7 #no unit
    BSFC = 0.48 #Brake Specific fuel consumption, no unit
    Payload = 40000 #in lb
    FS = 4 #Security factor, no unit
    Manufacturing_factor = 1.2 #no unit
    nb_septum = 2 #no unit,a voir si je modifie ca pour permettre de le modif
    Faf = 1.26 #account for the external attach fitting, no unit
    Seaming_factor = 1.06 #no unit
    nb_ball = nb_lobes*2 #6 #no unit
    Fpsq = 1.0 #account for the fabrication of the tail in lightweight structural material, no unit
    surface_CS = 0.2 #percentage of the CS on the tail surface, no unit
    Fact = 0.79 #no unit
    Finstall = 1.15 #account for the installation of actuators, no unit
    lec = 150 #account for the length of control cable for the engines, in ft
    Nb_prop = 4 #no unit
    Nbl = 3 #no unit
    Kp = 31.92 #no unit
    nt = 2 #number of fuel tank, no unit
    Integral_tank = 0 #percentage of fuel tank that are integral, no unit
    Npad = 3 #number of pads for ACLS, no unit
    Vsr = 4 #landing sink rate, in f/s
    Wcomputer = 250 #account for the weight of the computer, it's estimated, in lb
    Wavionics = 250 #account for the weight of instruments, in lb
    Wfs = 470 #account for the weight of fuel system, in lb
    Kelect = 33.73 #account for a coef that depend on the type of range flight, no unit
    np = 0.65 #propeller efficiency, no unit


    class Geometry(om.ExplicitComponent):
        def setup(self):
            self.add_input('V',val=1000000.0)
            self.add_input('FR',val=3.0)
            self.add_output('V23',val=1.0)
            self.add_output('de',val=3.0)
            self.add_output('lb',val=1.0)
            self.add_output('dc',val=1.0)
            self.add_output('w',val=1.0)
            self.add_output('ht',val=1.0)
            self.add_output('AR',val=1.0)
            self.add_output('Swet_body_lobed',val=1.0)
            self.add_output('SHT',val=1.0)
            self.add_output('SVT',val=1.0)
            self.add_output('Wg',val=1.0)
            self.add_output('Wh0',val=1.0)
            self.add_output('Wh1',val=1.0)
            self.add_output('total_fuel',val=1.0)
            self.add_output('Woe',val=1.0)
            self.add_output('Cd0',val=1.0)
            self.add_output('K',val=1.0)
            self.add_output('L', val=1.0)

        def setup_partials(self):
            # Finite difference all partials.
            self.declare_partials('*', '*', method='fd')

        def compute(self,inputs,outputs):
            V = inputs['V']
            FR = inputs['FR']

            outputs['V23'] = V23 = V**(2/3)
            outputs['de'] = de = ((6*V)/(FR*pi))**(1/3) #equivalent diameter of body of revolution
            outputs['lb'] = lb = FR*de
            outputs['dc'] = dc = (de/(-0.0178*nb_lobes**2+0.361*nb_lobes+0.575))
            outputs['w'] = w = (1+nb_lobes)*(dc/2)
            outputs['ht'] = ht = dc
            outputs['AR'] = AR = ((4*w**2)/(pi*lb*w))
            p = 1.6075
            Swet_body_ellipsoid = pi*((lb**p*w**p+lb**p*ht**p+w**p*ht**p)/3)**(1/p)
            if nb_lobes == 1:
                p_ellipsoid = 2*pi*(dc/2)
                p_lobes = 2*pi*(dc/2)
            if nb_lobes == 2:
                p_ellipsoid = 2.525*pi*(dc/2)
                p_lobes = (8*pi*(dc/2))/3
            if nb_lobes == 3:
                p_ellipsoid = 3.084*pi*(dc/2)
                p_lobes = (10*pi*(dc/2))/3
            if nb_lobes == 4:
                p_ellipsoid = 3.663*pi*(dc/2)
                p_lobes = (12*pi*(dc/2))/3
            if nb_lobes == 5:
                p_ellipsoid = 4.254*pi*(dc/2)
                p_lobes = (14*pi*(dc/2))/3
            outputs['Swet_body_lobed'] = Swet_body_lobed = (p_lobes/p_ellipsoid)*Swet_body_ellipsoid
            CHT = -0.0051*(1000000/V)+0.0717
            outputs['SHT'] = SHT = CHT*(V23*lb)/(0.38*lb)
            CVT = -0.0049*(1000000/V)+0.0641
            outputs['SVT'] = SVT = CVT*(V23*lb)/(0.38*lb)
            Cd0tab = []
            q = (1/2)*pho_cruise_alt*Vs**2
            Re = ((lb*Vs*pho_cruise_alt)/mu)
            try:
                Cf = (0.455/(log10(Re)**2.58))
            except ValueError as e:
                raise om.AnalysisError('Re out of range')
            FF3D_body = 1+(1.5/FR**1.5)+(7/FR**3)
            Cd0_body = ((FF3D_body*Cf*Swet_body_lobed)/V23)
            Cd0tab.append(Cd0_body)
            FF_tail = 1+1.2*tc+100*tc**4
            ctail = ((AR_tails*SHT*0.5)**0.5+(AR_tails*SVT*0.5)**0.5)/2
            Re_tail = (ctail*Vs*pho_cruise_alt)/mu
            cf_tail = 0.455/(log10(Re_tail)**2.58)
            Swet_tail = 2.2*(SHT+SVT)
            Cd0_tail = (FF_tail*cf_tail*Swet_tail)/V23
            Cd0tab.append(Cd0_tail)
            #Add every Cd0 for each part of the airship, check what we choose in the end
            Cd0cab = (0.108*Cd0_body*V23+7.7)/V23
            Cd0tab.append(Cd0cab)
            Cd0eng = NE*4.25/V23
            Cd0tab.append(Cd0eng)
            Cd0engcooling = NE*(2*10**(-6)*V+4.1)/V23
            Cd0tab.append(Cd0engcooling)
            Cd0engmount = (0.044*Cd0_body*V23+0.92)/V23
            Cd0tab.append(Cd0engmount)
            Cd0cables = (9.7*10**(-6)*V+10.22)/V23
            Cd0tab.append(Cd0cables)
            Cd0ACLS = 0.0002
            Cd0tab.append(Cd0ACLS)
            Cd0_int = (4.78*10**(-6)*V)/V23
            Cd0tab.append(Cd0_int)
            outputs['Cd0'] = Cd0 = sum(Cd0tab)
            outputs['K'] = -0.0145*(1/AR)**4 + 0.182*(1/AR)**3 - 0.514*(1/AR)**2 + 0.838*(1/AR) - 0.053
            if nb_lobes == 1:
                Nl = 2
            if nb_lobes == 2:
                Nl = 2.25
            if nb_lobes == 3:
                Nl = 2.4
            if nb_lobes == 4:
                Nl = 2.5
            if nb_lobes == 5:
                Nl = 2.54
            outputs['K'] = K = outputs['K']/Nl
            outputs['L'] = L = 0.0646*V*(pho_max_alt/pho_sea_level)
            Wlanding = L/BR
            outputs['Wh1'] = Wh1 = Wlanding - L
            A = (326*np)/(BSFC*sqrt(K*Cd0))
            B = q*V23*sqrt(Cd0/K)
            outputs['Wh0'] = Wh0 = B*tan((range/A)+atan(Wh1/B))
            fuel_burned = Wh0-Wh1
            fuelres = fuel_burned*0.05
            outputs['total_fuel'] = fuel_burned+fuelres
            Wzf = (L/BR)-fuelres
            outputs['Woe'] = Wzf - Payload
            outputs['Wg'] = Wlanding + fuel_burned

    class aerodynamic(om.ExplicitComponent):
        def setup(self):
            self.add_input('Wh0',val=1.0)
            self.add_input('V23',val=1.0)
            self.add_input('Cd0',val=1.0)
            self.add_input('K',val=1.0)
            self.add_output('qmax',val=1.0)
            self.add_output('Maxpower_perengine',val=1.0)
            self.add_output('Dp',val=1.0)
            self.add_output('np_check',val=1.0)

        def setup_partials(self):
            self.declare_partials('*','*',method='fd')

        def compute(self,inputs,outputs):
            Wh0= inputs['Wh0']
            V23= inputs['V23']
            Cd0= inputs['Cd0']
            K= inputs['K']
            qmax = outputs['qmax']
            Maxpower_perengine = outputs['Maxpower_perengine']
            Dp= outputs['Dp']
            np_check= outputs['np_check']


            L_aero = Wh0
            qmax[...] = pho_sea_level*(((Vs*1.1)**2)/2)
            Clmax_power = (L_aero/qmax/V23)
            Drag = (Cd0+K*Clmax_power**2)*qmax*V23
            Maxpower_perengine[...] = (Vs*1.1)*(Drag/np/NE/550)
            Cs = ((pho_sea_level*Vs**5)/Maxpower_perengine/550/10**2)**(1/5) #10 correspond a une rotation max pour le propeller speed and pho is in slug/ft3
            J = 0.156*Cs**2+0.241*Cs+0.138
            Dp[...] = ((Vs)/10)/J
            np_check[...] = 0.139*Cs**3-0.749*Cs**2+1.37*Cs+0.0115

    class second_weight_estimation(om.ExplicitComponent):
        def setup(self):
            self.add_input('qmax',val=1.0)
            self.add_input('ht',val=1.0)
            self.add_input('V',val=1000000.0)
            self.add_input('SHT',val=1.0)
            self.add_input('SVT',val=1.0)
            self.add_input('Swet_body_lobed',val=1.0)
            self.add_input('dc',val=1.0)
            self.add_input('lb',val=1.0)
            self.add_input('Maxpower_perengine',val=1.0)
            self.add_input('Dp',val=1.0)
            self.add_input('total_fuel',val=1.0)
            self.add_input('Woe',val=1.0)
            self.add_input('Wh0',val=1.0)
            self.add_input('Wh1',val=1.0)
            self.add_output('Wg2',val=1.0)

        def setup_partials(self):
            self.declare_partials('*','*',method='fd')

        def compute(self,inputs,outputs):
            qmax = inputs['qmax']
            ht = inputs['ht']
            V = inputs['V']
            SHT = inputs['SHT']
            SVT = inputs['SVT']
            Swet_body_lobed = inputs['Swet_body_lobed']
            dc = inputs['dc']
            lb = inputs['lb']
            Maxpower_perengine = inputs['Maxpower_perengine']
            Dp = inputs['Dp']
            total_fuel = inputs['total_fuel']
            Woe = inputs['Woe']
            Wh0 = inputs['Wh0']
            Wh1 = inputs['Wh1']
            Wg2 = outputs['Wg2']

            Weight = []
            Pint = 1.2*qmax+0.0635*ht
            hull_fabric_load = FS*12*(Pint/144)*(dc/2)
            hull_fabric_density = 0.0085*hull_fabric_load+1.365
            Whull = hull_fabric_density*Manufacturing_factor*Faf*Swet_body_lobed/(16*9) #to convert each unit in yd2 and then go into lb
            Septum_fabric_load = 1.5*hull_fabric_load
            Septum_fabric_density = 0.0085*Septum_fabric_load+1.365
            side_body_area = 0.75 #percentage of envelope side body area
            Wseptum = nb_septum*Seaming_factor*Septum_fabric_density*side_body_area*pi*ht*(lb/4)/(16*9) #we have a formula for the body side area after side body area
            Wenv = Wseptum+Whull
            Weight.append(Wenv)
            Vball = V*((1/(pho_max_alt/pho_sea_level))-1)
            Sball = nb_ball*pi*(nb_lobes*((Vball/pi)/6))**(2/3)
            Wball = 0.035*Sball
            Weight.append(Wball)
            Wssf = Fpsq*(SHT+SVT)*Faf*(1-surface_CS)
            Wcs = Fpsq*(SHT+SVT)*surface_CS
            Wtail = Wcs+Wssf
            Weight.append(Wtail)
            Wact = Fact*Finstall*surface_CS*(SHT+SVT)
            Weng = NE*Maxpower_perengine**(0.7956)*4.848
            Weight.append(Weng)
            Wengmt = 0.64*Weng
            Weight.append(Wengmt)
            Wec = 60.27*(lec*(NE/100))**(0.724)
            Weight.append(Wec)
            Wstart = 50.38*((Weng/1000))**(0.459)
            Weight.append(Wstart)
            Wprop = Nb_prop*Kp*Nbl**(0.391)*(Dp*(Maxpower_perengine/1000))**(0.782)
            Weight.append(Wprop)
            Wft = 2.49*(total_fuel/6)**(0.6)*nt**(0.2)*NE**(0.13)*(1/(1+Integral_tank))**(0.3) #total fuel is divided by 6 which the weight in lb per gallon for aviation gas
            Weight.append(Wft)
            Wpress_syst = 0.02*Woe
            Weight.append(Wpress_syst)
            Ppad = Pint
            ACLS_area_landing = 3*(0.23*Wh1*(Vsr/(Npad*Ppad))) #we multiply it by 3 because 3 pads but the configuration can change some stuff and then modify the calculations
            ACLS_area_takeoff = Wh0/Ppad
            Wacls = 1.6*max(ACLS_area_landing,ACLS_area_takeoff)
            Weight.append(Wacls)
            Wvms = Wcomputer+Wavionics+Wact
            Weight.append(Wvms)
            Welect = Kelect * (Wfs+Wcomputer+Wavionics)**(0.51)  #maybe don't have this equation
            Weight.append(Welect)
            Wmsys = 0.05*Woe
            Weight.append(Wmsys)
            Wuf = 0.01*total_fuel
            Weight.append(Wuf)
            Wmargin = 0.06*Woe
            Weight.append(Wmargin)
            Woe2 = sum(Weight)
            Wg2[...] = Woe2 + total_fuel + Payload

    class Optimization(om.Group):
        def setup(self):
            cycle = self.add_subsystem('cycle',om.Group(),
                                       promotes_inputs=['V','FR'],
                                       promotes_outputs=['Wg', 'Wg2', 'L', 'np_check'])
            cycle.add_subsystem('d1', Geometry(),
                                promotes_inputs=['V','FR'],
                                promotes_outputs=['Wg', 'L'])
            cycle.add_subsystem('d4', aerodynamic(),
                                promotes_outputs=['np_check'])
            cycle.add_subsystem('d5', second_weight_estimation(),
                                promotes_inputs=['V'],
                                promotes_outputs=['Wg2'])

            cycle.connect('d1.V23','d4.V23')
            cycle.connect('d1.lb','d5.lb')
            cycle.connect('d1.dc','d5.dc')
            cycle.connect('d1.ht','d5.ht')
            cycle.connect('d1.Swet_body_lobed','d5.Swet_body_lobed')
            cycle.connect('d1.SHT','d5.SHT')
            cycle.connect('d1.SVT','d5.SVT')
            cycle.connect('d1.Wh0','d4.Wh0')
            cycle.connect('d1.Wh0','d5.Wh0')
            cycle.connect('d1.Wh1','d5.Wh1')
            cycle.connect('d1.total_fuel','d5.total_fuel')
            cycle.connect('d1.Woe','d5.Woe')
            cycle.connect('d1.Cd0','d4.Cd0')
            cycle.connect('d1.K','d4.K')
            cycle.connect('d4.qmax','d5.qmax')
            cycle.connect('d4.Maxpower_perengine','d5.Maxpower_perengine')
            cycle.connect('d4.Dp','d5.Dp')

            # No feedback, not necessary
            # cycle.nonlinear_solver = om.NonlinearBlockGS()

            cycle.set_input_defaults('V',val=2)

            self.add_subsystem('obj', om.ExecComp('obj = Wg - Wg2'), promotes= ['obj','Wg','Wg2'])
            self.add_subsystem('c1', om.ExecComp('c1 = L/Wg'), promotes=['c1','L','Wg'])
            self.add_subsystem('c2', om.ExecComp('c2 = np_check'), promotes=['c2','np_check'])


    prob = om.Problem()
    prob.model = Optimization()

    prob.driver = om.ScipyOptimizeDriver()
    prob.driver.options['maxiter'] = 100
    prob.driver.options['tol'] = 1e-8

    prob.model.add_design_var('V', lower=0, ref=1.0E6)
    prob.model.add_design_var('FR', lower=0)
    prob.model.add_objective('obj', ref=1.0E6)
    prob.model.add_constraint('c1', lower=0, upper=1.0)
    prob.model.add_constraint('np_check', lower=0.95*np, upper=1.05*np)

    # Ask OpenMDAO to finite-difference across the model to compute the gradients for the optimizer
    prob.model.approx_totals()

    prob.setup()
    prob.set_val('V',2000000)
    prob.set_val('FR',3)

    prob.run_driver()

此时，我认为 OpenMDAO 问题的结构是正常的（事物似乎连接正确），但由于雷诺数问题，有关优化问题的某些内容无法正常工作。

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.