diff --git a/gpkitmodels/SP/SimPleAC/README.md b/gpkitmodels/SP/SimPleAC/README.md index 66c44493..c8b59e1d 100644 --- a/gpkitmodels/SP/SimPleAC/README.md +++ b/gpkitmodels/SP/SimPleAC/README.md @@ -290,3 +290,8 @@ We have tested a variety of potential objectives for the SimpPleAC model, some o \item $\frac{W_f}{T_{flight}}$: Product of the fuel weight and the inverse of the time of flight. \item $W_{f} + c \times T_{flight}$: A linear combination of fuel weight and time of flight. This can simulate recurring costs (fuel and labor), and yield interesting results. \end{itemize} + +Extension to mission and multimission design +------------ + +This model has been developed extensively in [Berk's Master's thesis](http://convex.mit.edu/publications/ozturk_masters_thesis.pdf) where it has been extended to [single-mission](https://github.com/convexengineering/gplibrary/blob/master/gpkitmodels/SP/SimPleAC/SimPleAC_mission.py) and [multi-mission](https://github.com/convexengineering/gplibrary/blob/master/gpkitmodels/SP/SimPleAC/SimPleAC_multimission.py) scenarios. diff --git a/gpkitmodels/SP/SimPleAC/SimPleAC_mission4.py b/gpkitmodels/SP/SimPleAC/SimPleAC_mission4.py new file mode 100644 index 00000000..574ee90e --- /dev/null +++ b/gpkitmodels/SP/SimPleAC/SimPleAC_mission4.py @@ -0,0 +1,290 @@ +from gpkit import Model, Variable, SignomialsEnabled, SignomialEquality, VarKey, units,Vectorize +from gpkit.constraints.bounded import Bounded +from gpkit.constraints.tight import Tight +from relaxed_constants import relaxed_constants, post_process + +import numpy as np +import matplotlib.pyplot as plt + +# Importing models +from atmosphere import Atmosphere + +# SimPleAC with mission design and flight segments, and lapse rate and BSFC model (3.4.2) + +class SimPleAC(Model): + def setup(self): + self.engine = Engine() + self.wing = Wing() + self.fuse = Fuselage() + self.components = [self.engine, self.wing, self.fuse] + + # Environmental constants + g = Variable("g", 9.81, "m/s^2", "gravitational acceleration") + rho_f = Variable("\\rho_f", 817, "kg/m^3", "density of fuel") + + # Free Variables + W = Variable("W", "N", "maximum takeoff weight") + W_f = Variable("W_f", "N", "maximum fuel weight") + V_f = Variable("V_f", "m^3", "maximum fuel volume") + V_f_avail = Variable("V_{f_{avail}}", "m^3", "fuel volume available") + + constraints = [] + + # Thrust and drag model + constraints += [self.fuse['C_{D_{fuse}}'] == self.fuse['(CDA0)'] / self.wing['S']] + + # Fuel volume model + with SignomialsEnabled(): + constraints += [V_f == W_f / g / rho_f, + V_f_avail <= self.wing['V_{f_{wing}}'] + self.fuse['V_{f_{fuse}}'], #[SP] + V_f_avail >= V_f] + + return constraints, self.components + + def dynamic(self,state): + return SimPleACP(self,state) + +class SimPleACP(Model): + def setup(self,aircraft,state): + self.aircraft = aircraft + self.engineP = aircraft.engine.dynamic(state) + self.wingP = aircraft.wing.dynamic(state) + self.Pmodels = [self.engineP, self.wingP] + + # Free variables + C_D = Variable("C_D", "-", "drag coefficient") + D = Variable("D", "N", "total drag force") + LoD = Variable('L/D','-','lift-to-drag ratio') + Re = Variable("Re", "-", "Reynolds number") + V = Variable("V", "m/s", "cruising speed") + + constraints = [] + + constraints += [self.engineP['T'] * V <= self.aircraft.engine['\\eta_{prop}'] * self.engineP['P_{shaft}'], + C_D >= self.aircraft['C_{D_{fuse}}'] + self.wingP['C_{D_{wpar}}'] + self.wingP['C_{D_{ind}}'], + D >= 0.5 * state['\\rho'] * self.aircraft['S'] * C_D * V ** 2, + Re == (state['\\rho'] / state['\\mu']) * V * (self.aircraft['S'] / self.aircraft['A']) ** 0.5, + self.wingP['C_f'] >= 0.074 / Re ** 0.2, + LoD == self.wingP['C_L'] / C_D] + + return constraints, self.Pmodels + +class Fuselage(Model): + def setup(self): + # Free Variables + CDA0 = Variable("(CDA0)", "m^2", "fuselage drag area") #0.035 originally + C_D_fuse = Variable('C_{D_{fuse}}','-','fuselage drag coefficient') + # Free variables (fixed for performance eval.) + V_f_fuse = Variable('V_{f_{fuse}}','m^3','fuel volume in the fuselage', fix = True) + + constraints = [] + constraints += [V_f_fuse == 10*units('m')*CDA0, + V_f_fuse >= 1*10**-5*units('m^3')] + + return constraints + +class Wing(Model): + def setup(self): + # Non-dimensional constants + C_Lmax = Variable("C_{L,max}", 1.6, "-", "lift coefficient at stall", pr=5.) + e = Variable("e", 0.92, "-", "Oswald efficiency factor", pr=3.) + k = Variable("k", 1.17, "-", "form factor", pr=10.) + N_ult = Variable("N_{ult}", 3.3, "-", "ultimate load factor", pr=15.) + S_wetratio = Variable("(\\frac{S}{S_{wet}})", 2.075, "-", "wetted area ratio", pr=3.) + tau = Variable("\\tau", 0.12, "-", "airfoil thickness to chord ratio", pr=10.) + + # Dimensional constants + W_w_coeff1 = Variable("W_{w_{coeff1}}", 2e-5, "1/m", + "wing weight coefficient 1", pr= 30.) #orig 12e-5 + W_w_coeff2 = Variable("W_{w_{coeff2}}", 60., "Pa", + "wing weight coefficient 2", pr=10.) + + # Free Variables (fixed for performance eval.) + A = Variable("A", "-", "aspect ratio",fix = True) + S = Variable("S", "m^2", "total wing area", fix = True) + W_w = Variable("W_w", "N", "wing weight")#, fix = True) + W_w_strc = Variable('W_{w_{strc}}','N','wing structural weight', fix = True) + W_w_surf = Variable('W_{w_{surf}}','N','wing skin weight', fix = True) + V_f_wing = Variable("V_{f_{wing}}",'m^3','fuel volume in the wing', fix = True) + + constraints = [] + + # Structural model + constraints += [W_w_surf >= W_w_coeff2 * S, + W_w >= W_w_surf + W_w_strc] + + # Wing fuel model + constraints += [V_f_wing**2 <= 0.0009*S**3/A*tau**2] # linear with b and tau, quadratic with chord + + return constraints + + def dynamic(self,state): + return WingP(self,state) + +class WingP(Model): + def setup(self,wing,state): + self.wing = wing + # Free Variables + C_f = Variable("C_f", "-", "skin friction coefficient") + C_D_ind = Variable('C_{D_{ind}}', '-', "wing induced drag") + C_D_wpar = Variable('C_{D_{wpar}}', '-', 'wing profile drag') + C_L = Variable("C_L", "-", "wing lift coefficient") + + constraints = [] + + # Drag model + constraints += [C_D_ind == C_L ** 2 / (np.pi * self.wing['A'] * self.wing['e']), + C_D_wpar == self.wing['k'] * C_f * self.wing["(\\frac{S}{S_{wet}})"]] + + return constraints + +class Engine(Model): + def setup(self): + # Dimensional constants + BSFC_ref = Variable("BSFC_{ref}", 0.32, "lbf/(hp*hr)", "reference brake specific fuel consumption") + eta_prop = Variable("\\eta_{prop}",0.8,'-',"propeller efficiency") + P_shaft_ref = Variable("P_{shaft,ref}",149,"kW","reference MSL maximum shaft power") + W_e_ref = Variable("W_{e,ref}",153, "lbf","reference engine weight") + h_ref = Variable("h_{ref}",15000,'ft','engine lapse reference altitude') + + # Free variables + P_shaft_max = Variable("P_{shaft,max}","kW","MSL maximum shaft power") + W_e = Variable("W_e","N","engine weight") + + constraints = [] + constraints += [(W_e/W_e_ref)**1.92 >= 0.00441 * (P_shaft_max/P_shaft_ref)**0.759 + + 1.44 * (P_shaft_max/P_shaft_ref)**2.90] + return constraints + + def dynamic(self,state): + return EngineP(self,state) + +class EngineP(Model): + def setup(self,engine,state): + self.engine = engine + # Dimensional constants + + # Free variables + BSFC = Variable("BSFC", "lbf/(hp*hr)", "brake specific fuel consumption") + P_shaft = Variable("P_{shaft}","kW","shaft power") + P_shaft_alt = Variable("P_{shaft,alt}","kW",'maximum shaft power at altitude') + Thrust = Variable("T","N","propeller thrust") + + L = Variable("L","-","power lapse percentage") + + constraints = [] + + with SignomialsEnabled(): + constraints += [P_shaft <= P_shaft_alt, + L == (0.937 * (state['h']/self.engine['h_{ref}'])**0.0922)**10, + SignomialEquality(1, L + P_shaft_alt / self.engine['P_{shaft,max}']), + (BSFC/self.engine['BSFC_{ref}'])**(0.1) >= 0.984*(P_shaft/P_shaft_alt)**-0.0346, + BSFC/self.engine['BSFC_{ref}'] >= 1., + ] + return constraints + + +class Mission(Model): + def setup(self,aircraft,Nsegments): + self.aircraft = aircraft + W_f_m = Variable('W_{f_m}','N','total mission fuel') + t_m = Variable('t_m','hr','total mission time') + + with Vectorize(Nsegments): + Wavg = Variable('W_{avg}','N','segment average weight') + Wstart = Variable('W_{start}', 'N', 'weight at the beginning of flight segment') + Wend = Variable('W_{end}', 'N', 'weight at the end of flight segment') + h = Variable('h','m','final segment flight altitude') + havg = Variable('h_{avg}','m','average segment flight altitude') + dhdt = Variable('\\frac{dh}{dt}','m/hr','climb rate') + W_f_s = Variable('W_{f_s}','N', 'segment fuel burn') + t_s = Variable('t_s','hr','time spent in flight segment') + R_s = Variable('R_s','km','range flown in segment') + state = Atmosphere() + self.aircraftP = self.aircraft.dynamic(state) + + # Mission variables + hcruise = Variable('h_{cruise_m}', 'm', 'minimum cruise altitude') + Range = Variable("Range_m", "km", "aircraft range") + W_p = Variable("W_{p_m}", "N", "payload weight", pr=20.) + V_min = Variable("V_{min_m}", "m/s", "takeoff speed", pr=20.) + TOfac = Variable('T/O factor_m', '-','takeoff thrust factor') + cost_index = Variable("C_m", '1/hr','hourly cost index') + + constraints = [] + + # Setting up the mission + with SignomialsEnabled(): + constraints += [havg == state['h'], # Linking states + h[1:Nsegments-1] >= hcruise, # Adding minimum cruise altitude + + # Weights at beginning and end of mission + Wstart[0] >= W_p + self.aircraft.wing['W_w'] + self.aircraft.engine['W_e'] + W_f_m, + Wend[Nsegments-1] >= W_p + self.aircraft.wing['W_w'] + self.aircraft.engine['W_e'], + + # Lift, and linking segment start and end weights + Wavg <= 0.5 * state['\\rho'] * self.aircraft['S'] * self.aircraftP.wingP['C_L'] * self.aircraftP['V'] ** 2, + Wstart >= Wend + W_f_s, # Making sure fuel gets burnt! + Wstart[1:Nsegments] == Wend[:Nsegments-1], + Wavg == Wstart ** 0.5 * Wend ** 0.5, + + # Altitude changes + h[0] == t_s[0]*dhdt[0], # Starting altitude + dhdt >= 1.*units('m/hr'), + havg[0] == 0.5*h[0], + havg[1:Nsegments] == (h[1:Nsegments]*h[0:Nsegments-1])**(0.5), + SignomialEquality(h[1:Nsegments],h[:Nsegments-1] + t_s[1:Nsegments]*dhdt[1:Nsegments]), + + # Thrust and fuel burn + W_f_s >= self.aircraftP.engineP['BSFC'] * self.aircraftP.engineP['P_{shaft}'] * t_s, + self.aircraftP.engineP['T'] * self.aircraftP['V'] >= self.aircraftP['D'] * self.aircraftP['V'] + Wavg * dhdt, + + # Max MSL thrust at least 2*climb thrust + self.aircraft.engine['P_{shaft,max}'] >= TOfac*self.aircraftP.engineP['P_{shaft}'][0], + # Flight time + t_s == R_s/self.aircraftP['V'], + + # Aggregating segment variables + self.aircraft['W_f'] >= W_f_m, + R_s == Range/Nsegments, # Dividing into equal range segments + W_f_m >= sum(W_f_s), + t_m >= sum(t_s) + ] + + # Maximum takeoff weight + constraints += [self.aircraft['W'] >= W_p + self.aircraft.wing['W_w'] + self.aircraft['W_f'] + self.aircraft.engine['W_e']] + + # Stall constraint + constraints += [self.aircraft['W'] <= 0.5 * state['\\rho'] * + self.aircraft['S'] * self.aircraft['C_{L,max}'] * V_min ** 2] + + # Wing weight model + constraints += [self.aircraft.wing['W_{w_{strc}}']**2. >= + self.aircraft.wing['W_{w_{coeff1}}']**2. / self.aircraft.wing['\\tau']**2. * + (self.aircraft.wing['N_{ult}']**2. * self.aircraft.wing['A'] ** 3. * + ((W_p+self.aircraft.fuse['V_{f_{fuse}}']*self.aircraft['g']*self.aircraft['\\rho_f']) * + self.aircraft['W'] * self.aircraft.wing['S']))] + + # Upper bounding variables + constraints += [t_m <= 100000*units('hr'), + W_f_m <= 1e10*units('N')] + + return constraints, state, self.aircraft, self.aircraftP + +if __name__ == "__main__": + # Most basic way to execute the model + m = Mission(SimPleAC(),4) + m.substitutions.update({ + 'h_{cruise_m}' :5000*units('m'), + 'Range_m' :3000*units('km'), + 'W_{p_m}' :6250*units('N'), + 'C_m' :120*units('1/hr'), + 'V_{min_m}' :25*units('m/s'), + 'T/O factor_m' :2, + }) + m.cost = m['W_{f_m}']*units('1/N') + m['C_m']*m['t_m'] + #m = Model(m.cost, Bounded(m)) + #m_relax = relaxed_constants(m,None,None) + sol = m.localsolve(verbosity = 4) + #post_process(sol) + #print sol.table() diff --git a/gpkitmodels/SP/SimPleAC/SimPleAC_multimission.py b/gpkitmodels/SP/SimPleAC/SimPleAC_multimission.py new file mode 100644 index 00000000..9da231db --- /dev/null +++ b/gpkitmodels/SP/SimPleAC/SimPleAC_multimission.py @@ -0,0 +1,68 @@ +from gpkit import Model, Variable, SignomialsEnabled, SignomialEquality, VarKey, units,Vectorize +from gpkit.constraints.bounded import Bounded +from relaxed_constants import relaxed_constants, post_process + +import numpy as np +import matplotlib.pyplot as plt + +from SimPleAC_mission4 import Mission, SimPleAC +from atmosphere import Atmosphere + +# SimPleAC with multimission design (3.5) + +class Multimission(Model): + def setup(self,aircraft,Nmissions,Nsegments): + self.aircraft = aircraft + self.missions = [] + for i in range(0,Nmissions): + self.missions.append(Mission(self.aircraft,Nsegments)) + + # Multimission objective variables + W_f_mm = Variable('W_{f_{mm}}','N','multimission fuel weight') + + with Vectorize(Nmissions): + # Mission variables + hcruise = Variable('h_{cruise_{mm}}', 'm', 'minimum cruise altitude') + Range = Variable("Range_{mm}", "km", "aircraft range") + W_p = Variable("W_{p_{mm}}", "N", "payload weight", pr=20.) + V_min = Variable("V_{min_{mm}}", 25, "m/s", "takeoff speed", pr=20.) + cost_index = Variable("C_{mm}", '1/hr','hourly cost index') + TOfac = Variable('T/O factor_{mm}', 2.,'-','takeoff thrust factor') + + constraints = [] + + # Setting up the missions + for i in range(0,Nmissions): + constraints += [ + self.missions[i]['h_{cruise_m}'] == hcruise[i], + self.missions[i]['Range_m'] == Range[i], + self.missions[i]['W_{p_m}'] == W_p[i], + self.missions[i]['V_{min_m}'] == V_min[i], + self.missions[i]['C_m'] == cost_index[i], + self.missions[i]['T/O factor_m'] == TOfac[i], + # Upper bounding relevant variables + W_f_mm <= 1e11*units('N'), + ] + + # Multimission constraints + constraints += [W_f_mm >= sum(self.missions[i]['W_{f_m}'] for i in range(0,Nmissions))] + + + + return constraints, self.aircraft, self.missions + +if __name__ == "__main__": + Nmissions = 2 + Nsegments = 4 + aircraft = SimPleAC() + m = Multimission(aircraft,Nmissions,Nsegments) + m.substitutions.update({ + 'h_{cruise_{mm}}':[5000*units('m'), 5000*units('m')], + 'Range_{mm}' :[3000*units('km'), 2000*units('km')], + 'W_{p_{mm}}' :[6250*units('N'), 8000*units('N')], + 'C_{mm}' :[120*units('1/hr'), 360*units('1/hr')], + }) + #m.cost = m['W_{f_{mm}}']*units('1/N') + sum(m.missions[i]['C_m']*m.missions[i]['t_m'] for i in range(0,Nmissions)) + m.cost = (m.missions[0]['W_{f_m}']*units('1/N') + m.missions[1]['C_m']*m.missions[1]['t_m']) + sol = m.localsolve(verbosity = 4) +