在 OpenMDAO 中创建约束

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

欢迎大家阅读这篇文章,

所以我尝试使用模型的 2 个变量创建约束,但似乎我不知道如何设置它。

我试图确保 np_check 包含在 np 的 95% 和 105% 之间,其中 np 是一个设计变量,但它不知道我是否必须在优化类中或外部进行设置。

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=2000000.0)
        self.add_input('FR',val=3.0)
        self.add_input('np',val = 0.65)
        self.add_output('V23',val=15874.0)
        self.add_output('de',val=108.4)
        self.add_output('lb',val=325.2)
        self.add_output('dc',val=72.3)
        self.add_output('w',val=144.5)
        self.add_output('ht',val=72.3)
        self.add_output('AR',val=0.566)
        self.add_output('Swet_body_lobed',val=100101.0)
        self.add_output('SHT',val=2889.0)
        self.add_output('SVT',val=2575.0)
        self.add_output('Wg',val=170096.0)
        self.add_output('Wh0',val=68545.0)
        self.add_output('Wh1',val=43522.0)
        self.add_output('total_fuel',val=26274.0)
        self.add_output('Woe',val=103822.0)
        self.add_output('Cd0',val=0.03133)
        self.add_output('K',val=0.286)
        self.add_output('L', val=101551.0)
        self.add_output('qmax',val=10.26)
        self.add_output('Maxpower_perengine',val=867.0)
        self.add_output('Dp',val=21.2)
        self.add_output('np_check',val=0.667)

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

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

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

        qmax = outputs['qmax']
        Maxpower_perengine = outputs['Maxpower_perengine']
        Dp= outputs['Dp']
        np_check= outputs['np_check']

        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
        print('K',K)
        print('Cd0',Cd0)
        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
        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=10.26)
        self.add_input('ht',val=72.3)
        self.add_input('V',val=2000000.0)
        self.add_input('SHT',val=2889.0)
        self.add_input('SVT',val=2575.0)
        self.add_input('Swet_body_lobed',val=100101.0)
        self.add_input('dc',val=72.3)
        self.add_input('lb',val=325.2)
        self.add_input('Maxpower_perengine',val=867.0)
        self.add_input('Dp',val=21.2)
        self.add_input('total_fuel',val=26274.0)
        self.add_input('Woe',val=103822.0)
        self.add_input('Wh0',val=68545.0)
        self.add_input('Wh1',val=43522.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','np'],
                                    promotes_outputs=['Wg', 'Wg2', 'L', 'np_check'])
        cycle.add_subsystem('d1', Geometry(),
                            promotes_inputs=['V','FR','np'],
                            promotes_outputs=['Wg', 'L','np_check'])
        cycle.add_subsystem('d5', second_weight_estimation(),
                            promotes_inputs=['V'],
                            promotes_outputs=['Wg2'])

        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','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.qmax','d5.qmax')
        cycle.connect('d1.Maxpower_perengine','d5.Maxpower_perengine')
        cycle.connect('d1.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 = (Wg2 - Wg)**2'), promotes= ['obj','Wg','Wg2'])
        self.add_subsystem('c1', om.ExecComp('c1 = L/Wg'), promotes=['c1','L','Wg'])
        self.add_subsystem('lower_np', om.ExecComp('lower_np=0.95*np'), promotes=['lower_np','np'])
        self.add_subsystem('upper_np', om.ExecComp('upper_np=0.95*np'), promotes=['upper_np','np'])
        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=1000)
prob.model.add_design_var('FR', lower=1.0, upper=8.0)
prob.model.add_design_var('np', lower = 0.1)
prob.model.add_objective('obj',ref=1.0E3)
prob.model.add_constraint('c1', lower=0.1, upper=1.0)
prob.model.add_constraint('c2', lower='lower_np', upper='upper_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()
print('minimum found at')
print(prob.get_val('V'))
print(prob.get_val('FR'))
print(prob.get_val('Wg2'))
print(prob.get_val('np_check'))

print('minumum objective')
print(prob.get_val('obj')[0])

感谢您的阅读

optimization constraints openmdao
1个回答
0
投票

OpenMDAO 需要约束上的文字界限,因此您需要计算另一个输出:

np_check_con = np_check / np

因此,这将 np_check 计算为 np 的百分比。

然后添加约束

self.add_constraint('np_check_con', lower=0.95, upper=1.05)
© www.soinside.com 2019 - 2024. All rights reserved.