Welcome to everyone reading this post,
So I am trying to create a constraint using 2 variables of my model but it seems I don’t know how to setup this.
I am trying to make sure that np_check is contained between 95% of np and 105% of np where np is a design variable but it don’t know if I have to setup this in my Optimization class or outside.
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])
Thanks for reading