/
ssbj_idf.py
119 lines (102 loc) · 3.93 KB
/
ssbj_idf.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
"""
SSBJ test case - http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19980234657.pdf
Python implementation and OpenMDAO integration developed by
Sylvain Dubreuil and Remi Lafage of ONERA, the French Aerospace Lab.
"""
from __future__ import print_function
from sys import argv
import re
import numpy as np
import matplotlib.pylab as plt
from openmdao.api import Problem
from openmdao.api import SqliteRecorder, ScipyOptimizeDriver #, pyOptSparseDriver
from ssbj_idf_mda import SSBJ_IDF_MDA
from ssbj_mda import init_ssbj_mda
# pylint: disable=C0103
# Optimization problem
scalers = init_ssbj_mda()
print(scalers)
prob = Problem()
prob.model = SSBJ_IDF_MDA(scalers)
# Optimizer options
prob.driver = ScipyOptimizeDriver()
# prob.driver = pyOptSparseDriver()
optimizer = 'SLSQP'
prob.driver.options['optimizer'] = optimizer
#prob.driver.options['tol'] = 1e-3
#prob.driver.options['debug_print'] = ['desvars','ln_cons','nl_cons','objs']
#Design variables
prob.model.add_design_var('z', lower=np.array([0.2, 0.666,0.875,0.45,0.72,0.5]),
upper=np.array([1.8,1.333,1.125,1.45,1.27,1.5]))
prob.model.add_design_var('x_str', lower=np.array([0.4,0.75]), upper=np.array([1.6,1.25]))
prob.model.add_design_var('x_aer', lower=0.75, upper=1.25)
prob.model.add_design_var('x_pro', lower=0.18, upper=1.81)
prob.model.add_design_var('Theta')
prob.model.add_design_var('L')
prob.model.add_design_var('WE')
prob.model.add_design_var('WT')
prob.model.add_design_var('ESF')
prob.model.add_design_var('D')
# Objective function
prob.model.add_objective('obj')
#Constraints
prob.model.add_constraint('con_dt', upper=0.0)
prob.model.add_constraint('con_Theta_up', upper=0.0)
prob.model.add_constraint('con_Theta_low', upper=0.0)
prob.model.add_constraint('con_sigma1', upper=0.0)
prob.model.add_constraint('con_sigma2', upper=0.0)
prob.model.add_constraint('con_sigma3', upper=0.0)
prob.model.add_constraint('con_sigma4', upper=0.0)
prob.model.add_constraint('con_sigma5', upper=0.0)
prob.model.add_constraint('con_dpdx', upper=0.0)
prob.model.add_constraint('con_esf', upper=0.0)
prob.model.add_constraint('con_temp', upper=0.0)
#Coupling constraints
#Threshold for the coupling (constraints define as (x_in-x_out)**2<epsilon)
epsilon=1e-10
prob.model.add_constraint('con_str_aer_wt',upper=epsilon)
prob.model.add_constraint('con_str_aer_theta',upper=epsilon)
prob.model.add_constraint('con_aer_str_l',upper=epsilon)
prob.model.add_constraint('con_aer_pro_d',upper=epsilon)
prob.model.add_constraint('con_pro_aer_esf',upper=epsilon)
prob.model.add_constraint('con_pro_str_we',upper=epsilon)
#Recorder
db_name = 'ssbj_idf.sqlite'
if "--plot" in argv:
recorder = SqliteRecorder(db_name)
# recorder2 = WebRecorder("", case_name="SSBJ IDF")
prob.driver.recording_options['record_desvars'] = True
prob.driver.recording_options['record_objectives'] = True
prob.driver.recording_options['record_constraints'] = True
prob.driver.add_recorder(recorder)
# prob.driver.add_recorder(recorder2)
#Run optimization
prob.setup(mode='fwd')
prob.run_driver()
#prob.run_model()
#prob.check_partials()
#prob.cleanup()
print('Z_opt=', prob['z']*scalers['z'])
print('X_str_opt=', prob['x_str']*scalers['x_str'])
print('X_aer_opt=', prob['x_aer'])
print('X_pro_opt=', prob['x_pro']*scalers['x_pro'])
print('R_opt=', -prob['obj']*scalers['R'])
if "--plot" in argv:
from openmdao.recorders.case_reader import CaseReader
plt.figure()
cr = CaseReader(db_name)
print('Number of driver cases recorded =', cr.driver_cases.num_cases )
case_keys = cr.driver_cases.list_cases()
r = []
for case_key in case_keys:
r.append(-cr.driver_cases.get_case(case_key).objectives['obj']*scalers['R'])
plt.plot(r)
plt.xlabel('Iteration')
plt.ylabel('Range (Nm)')
plt.show()
# Check R =~ 3964Nm
R = float(-prob['obj']*scalers['R'])
assert(R > 3950.)
assert(R < 3970.)
# from openmdao.devtools.problem_viewer.problem_viewer import view_model
# view_model(prob)