Simon Boudoux
Reputation: 21
My design variables have no impact on the constraints and objective function in OpenMDAO
So I am working on an optimization problem with OpenMDAO. The subject of the optimization problem is to find the best Volume regarding some mission requirements that you chose and then make the Volume converge using 2 ways to compute the weight and compute the delta. So in the end my problem is after doing the code, I run it and then I got an error message saying that : "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." and "DerivativesWarning:Design variables [('V', inds=[0]), ('FR', inds=[0])] have no impact on the constraints or objective."
So in my code there is 3 subclass that compute the 2 weight and 1 about aerodynamics values. Then in my optimization class I put my 3 subsystem, I connect my inputs and outputs and then create my objective function and constraints and finally I make the problem work by implementing the optimization and then I try to run it and have my error message So here is the code (it has a lot of formulas but I think that is not really important in my subject I think ):
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'))
The total error message :
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.]
Thanks in advance for looking at it and hope that we can find the solution together
Upvotes: 0
Views: 124
Answers (1)
Rob Falck
Reputation: 2704
There are a few issues in the posted code.
When you write
np_check= outputs['np_check']
...
np_check = 5
Python does not write 5 into outputs['np_check'] as you might expect, but instead redefines it as a new integer. As a result, outputs['np_check'] maintains its default value (it's never overridden in the compute method).
To fix this, you can either assign to outputs['np_check'] directly, or force numpy to insert 5 into the array value of np_check using:
np_check[...] = 5
This explicitly tells numpy that you are shoving 5 into all elements of the np_check array (it only contains a single element in your instance).
After this, you'll still get an warning that the constraints and objectives cannot be impacted by your design variables. The best thing to do in this situation is to open up the n-squared diagram that OpenMDAO automatically creates for you in reports/n2.html of the directory where you executed the script.
In this diagram, variables are listed on the left side. Any variable colored orange is not connected to an output elsewhere in the model. For L and FR this makes sense because they are ultimately design variables provided by the optimizer.
But the diagram indicates that L, Wg, Wg2, and np_check are not connected to outputs. OpenMDAO is assuming that these inputs just hold whatever default value was assigned to them (typically 1.0), but the optimizer is not changing them.
These values appear to be outputs within the cycle group, but they need to be promoted, first out of their computing components, and then out of the cycle group itself.
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'])
FYI, I found this by clicking the "circled i" tool in the menubar on the left of the N2 diagram will allow you to see the "promoted" name of variables in the model. Variables which share the same promoted name are implicitly connected in OpenMDAO.
Another thing I did was to remove the nonlinear solver assignment from the cycle group. This group has no feedback, so the default nonlinear solver is preferred.
Finally, one of the constraints and the objective seemed to have values on the order of 1E6, so I scaled them by setting the ref value of those responses to 1.0E6. This will cause the optimizer to see the value of 1.0E6 as 1.0. There's a corresponding ref0 argument to add_objective and add_constraint that provides the value which the optimizer sees as zero.
With these changes, I get to a point where the optimizer is forcing the design into a region where log10(Re) is invalid. You can try raising an OpenMDAO AnalysisError when this occurs, which will force some optimizers to walk back to a valid region, but in this case it seems to insist on driving to that region.
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()
At this point I think the structure of the OpenMDAO problem is OK (things appear to be connected properly), but something about the optimization problem is not working correctly because of the Reynolds number issue.
Upvotes: 0
Related Questions
- Constraints or objectives cannot be impacted by the design variables of the problem
- How do I add design variables and constraints from AeroProblem and DVConstraints to OpenMDAO?
- openmdao compute gradient of constraints but not objective
- OpenMDAO: COBYLA and design variables boundaries
- openMDAO: Does the use of ExecComp maximum() interfere with constraints not being affected by design variables?
- Understanding the reasons behind Openmdao design
- OpenMDAO: Is there an optimum magnitude of each design variable or objective function to search design space?
- Discrete Independent/Design Variables ? openmdao 2.2.0
- OpenMDAO 1.x: recording desvars, constraints and objective
- OpenMDAO v0.13: Gradients of constraints and objectives