diff --git a/solar/npod_trade.py b/solar/npod_trade.py index 7ced97d..79539ca 100644 --- a/solar/npod_trade.py +++ b/solar/npod_trade.py @@ -1,24 +1,58 @@ " number of pods trade study " -from solar import Mission, Aircraft +# from solar import Mission, Aircraft import matplotlib.pyplot as plt import pandas as pd import numpy as np import sys +from relaxed_constants import relaxed_constants, post_process -def pods(N=[1, 3, 5, 7, 9, 0]): +def pods(N=[1, 3, 5, 7, 9, 0], Nplot=5): "trade number of pods" SP = True data = {} for i in N: + from solar import Mission, Aircraft Vehicle = Aircraft(Npod=i, sp=SP) M = Mission(Vehicle, latitude=[20]) M.cost = M[M.aircraft.Wtotal] - sol = M.localsolve("mosek") - data[i] = sol("Wtotal").magnitude + try: + sol = M.localsolve("mosek") + data[i] = sol("Wtotal").magnitude + except RuntimeWarning: + V2 = Aircraft(Npod=i, sp=SP) + M2 = Mission(V2, latitude=[20]) + M2.cost = M2[M2.aircraft.Wtotal] + feas = relaxed_constants(M2) + sol2 = feas.localsolve("mosek") + vks = post_process(sol2) + data[i] = np.NaN if vks else sol2("Wtotal").magnitude + M, sol = M2, sol2 + + if Nplot == i: + plot_shear(M, sol) df = pd.DataFrame(data, index=[0]) return df +def plot_shear(model, result): + " plot shear and moment diagrams " + + S = result(model.mission[1].winggust.S) + m = result(model.mission[1].winggust.M) + fig, ax = plt.subplots(2) + ax[0].plot(range(20), S) + ax[1].plot(range(20), m) + ax[0].grid(); ax[1].grid() + fig.savefig("shearandmoment.pdf") + + S = result(model.mission[1].wingg.S) + m = result(model.mission[1].wingg.M) + fig, ax = plt.subplots(2) + ax[0].plot(range(20), S) + ax[1].plot(range(20), m) + ax[0].grid(); ax[1].grid() + fig.savefig("shearandmoment2.pdf") + def plot_pods(df): " plot pod trade " @@ -44,7 +78,7 @@ def plot_pods(df): def test(): " for unit testing " - pods(N=[0]) + pods(Nplot=100) if __name__ == "__main__": if len(sys.argv) > 1: @@ -52,7 +86,7 @@ def test(): else: path = "" - GENERATE = False + GENERATE = True if GENERATE: DF = pods() diff --git a/solar/relaxed_constants.py b/solar/relaxed_constants.py new file mode 100644 index 0000000..b2f810d --- /dev/null +++ b/solar/relaxed_constants.py @@ -0,0 +1,51 @@ +from gpkit.constraints.relax import ConstantsRelaxed +from gpkit.constraints.bounded import Bounded +from gpkit import Model + +""" +Methods to precondition an SP so that it solves with a relaxed constants algorithim +and postcondition an SP to ensure all relax values are 1 +""" + +def relaxed_constants(model, include_only=None, exclude=None): + """ + Method to precondition an SP so it solves with a relaxed constants + algorithim + + ARGUMENTS + --------- + model: the model to solve with relaxed constants + + RETURNS + ------- + feas: the input model but with relaxed constants and a new objective + """ + + if model.substitutions: + constsrelaxed = ConstantsRelaxed(Bounded(model)) + feas = Model(constsrelaxed.relaxvars.prod()**20 * model.cost, + constsrelaxed) + # NOTE: It hasn't yet been seen but might be possible that + # the model.cost component above could cause infeasibility + else: + feas = Model(model.cost, model) + + return feas + +def post_process(sol): + """ + Model to print relevant info for a solved model with relaxed constants + + ARGUMENTS + -------- + sol: the solution to the solved model + """ + + bdvars = [d for d in sol.program.varkeys if "Relax" in str(d) + and "before" not in str(d) and sol(d).magnitude >= 1.001] + if bdvars: + print "GP iteration has relaxed constants" + print sol.program.result.table(varkeys) + + return bdvars + diff --git a/solar/sens_chart.py b/solar/sens_chart.py index 3d75b24..5f4049d 100644 --- a/solar/sens_chart.py +++ b/solar/sens_chart.py @@ -108,10 +108,10 @@ def test(): else: path = "" - V = Aircraft(Npod=3, sp=True) - M = Mission(V, latitude=[15]) + V = Aircraft(Npod=0, sp=False) + M = Mission(V, latitude=[20]) M.cost = M[M.aircraft.Wtotal] - sol = M.localsolve("mosek") + sol = M.solve("mosek") vns = {M.aircraft.Wpay: "$W_{\\mathrm{pay}}$", M.aircraft.battery.etacharge: "$\\eta_{\\mathrm{charge}}$", @@ -121,6 +121,6 @@ def test(): "Nmax": "$N_{\\mathrm{max}}$", "e": "$e$", "etaprop": "$\\eta_{\\mathrm{prop}}$"} - sd = get_highestsens(M, sol, N=15) + sd = get_highestsens(M, sol, N=10) f, a = plot_chart(sd) f.savefig(path + "sensbar.pdf", bbox_inches="tight") diff --git a/solar/solar.py b/solar/solar.py index 71f2845..608fcb9 100644 --- a/solar/solar.py +++ b/solar/solar.py @@ -5,16 +5,16 @@ from os import sep from numpy import hstack import pandas as pd -import numpy as np import gassolar.environment from ad.admath import exp from gassolar.environment.solar_irradiance import get_Eirr, twi_fits from gassolar.environment.wind_speeds import get_month -from gpkit import Model, parse_variables, Vectorize, Variable +from gpkit import Model, parse_variables, Vectorize, SignomialsEnabled from gpkit.tests.helpers import StdoutCaptured from gpkitmodels.GP.aircraft.wing.wing import Wing as WingGP from gpkitmodels.SP.aircraft.wing.wing import Wing as WingSP -from gpkitmodels.GP.aircraft.wing.boxspar import BoxSpar +from gpkitmodels.GP.aircraft.wing.boxspar import BoxSpar as BoxSparGP +from gpkitmodels.SP.aircraft.wing.boxspar import BoxSpar as BoxSparSP from gpkitmodels.GP.aircraft.wing.wing_skin import WingSecondStruct from gpkitmodels.GP.aircraft.tail.empennage import Empennage from gpkitmodels.GP.aircraft.tail.horizontal_tail import HorizontalTail @@ -28,6 +28,7 @@ from gpkitmodels.GP.aircraft.motor.motor import Motor from gpkitmodels import g from gpfit.fit_constraintset import FitCS as FCS +from relaxed_constants import relaxed_constants, post_process path = dirname(gassolar.environment.__file__) @@ -122,7 +123,7 @@ def setup(self, static, state, onDesign = False): self.wing.substitutions[e] = 0.95 self.wing.substitutions[self.wing.CLstall] = 4 - + self.wing.substitutions[e] = 0.95 dvars = [cdht*Sh/Sw, cdvt*Sv/Sw, cftb*Stb/Sw] @@ -163,8 +164,8 @@ class Aircraft(Model): minvttau 0.09 [-] minimum vertical tail tau ratio minhttau 0.06 [-] minimum horizontal tail tau ratio maxtau 0.144 [-] maximum wing tau ratio - - SKIP VERIFICATION + + SKIP VERIFICATION Upper Unbounded --------------- @@ -210,27 +211,27 @@ def setup(self, Npod=0, sp=False): foamhd.substitutions.update({foamhd.rho: 0.03}) materials = [cfrpud, cfrpfabric, foamhd] - HorizontalTail.sparModel = BoxSpar + HorizontalTail.sparModel = BoxSparGP HorizontalTail.fillModel = None HorizontalTail.skinModel = WingSecondStruct - VerticalTail.sparModel = BoxSpar + VerticalTail.sparModel = BoxSparGP VerticalTail.fillModel = None VerticalTail.skinModel = WingSecondStruct - TailBoom.__bases__ = (BoxSpar,) + TailBoom.__bases__ = (BoxSparGP,) TailBoom.secondaryWeight = True self.emp = Empennage(N=5) self.solarcells = SolarCells() self.battery = Battery() if sp: - WingSP.sparModel = BoxSpar + WingSP.sparModel = BoxSparSP WingSP.fillModel = None WingSP.skinModel = WingSecondStruct - self.wing = WingSP(N=10) + self.wing = WingSP(N=20) else: - WingGP.sparModel = BoxSpar + WingGP.sparModel = BoxSparGP WingGP.fillModel = None WingGP.skinModel = WingSecondStruct - self.wing = WingGP(N=10) + self.wing = WingGP(N=20) self.motor = Motor() Propeller.flight_model = ActuatorProp self.propeller = Propeller() @@ -462,23 +463,20 @@ def setup(self, aircraft, latitude=35, day=355): self.aircraft = aircraft self.fs = FlightState(latitude=latitude, day=day) - self.aircraftPerf = self.aircraft.flight_model(aircraft, self.fs, True) + self.aircraftPerf = self.aircraft.flight_model(aircraft, self.fs, False) self.slf = SteadyLevelFlight(self.fs, self.aircraft, self.aircraftPerf) - if aircraft.Npod is not 0: - if aircraft.Npod is not 1: - assert self.aircraft.sp - loadsp = self.aircraft.sp - else: - loadsp = False + if aircraft.Npod is not 0 and aircraft.Npod is not 1: + assert self.aircraft.sp + loadsp = self.aircraft.sp else: loadsp = False self.wingg = self.aircraft.wing.spar.loading( - self.aircraft.wing, self.fs, sp=loadsp) + self.aircraft.wing, self.fs, out=loadsp) self.winggust = self.aircraft.wing.spar.gustloading( - self.aircraft.wing, self.fs, sp=loadsp) + self.aircraft.wing, self.fs, out=loadsp) self.htailg = self.aircraft.emp.htail.spar.loading( self.aircraft.emp.htail, self.fs) self.vtailg = self.aircraft.emp.vtail.spar.loading( @@ -524,17 +522,38 @@ def setup(self, aircraft, latitude=35, day=355): self.vtailg.W == qne*Sv*CLvmax, ] - if self.aircraft.Npod is not 0: - if self.aircraft.Npod is not 1: - z0 = Variable("z0", 1e-10, "N", "placeholder zero value") - Nwing, Npod = self.aircraft.wing.N, self.aircraft.Npod - ypod = Nwing/((Npod-1)/2 + 1) - y0len = ypod-1 - weight = self.aircraft.battery.W/Npod*self.wingg.Nsafety - sout = np.hstack([[z0]*y0len + [weight]]*(Npod/2)) - sout = list(sout) + [z0]*(Nwing - len(sout) - 1) - constraints.extend([sout == self.wingg.Sout, - sout == self.winggust.Sout]) + if self.aircraft.Npod is not 0 and self.aircraft.Npod is not 1: + Nwing, Npod = self.aircraft.wing.N, self.aircraft.Npod + ypod = Nwing/((Npod-1)/2 + 1) + ypods = [ypod*n for n in range(1, (Npod-1)/2+1)] + Sgust, Mgust = self.winggust.S, self.winggust.M + qgust, Sg, Mg = self.winggust.q, self.wingg.S, self.wingg.M + qg = self.wingg.q + deta = self.aircraft.wing.planform.deta + b = self.aircraft.wing.planform.b + weight = self.aircraft.battery.W/Npod*self.wingg.N + for i in range(Nwing-1): + if i in ypods: + with SignomialsEnabled(): + constraints.extend([ + Sgust[i] >= (Sgust[i+1] + 0.5*deta[i]*(b/2) + * (qgust[i] + qgust[i+1]) - weight), + Sg[i] >= (Sg[i+1] + 0.5*deta[i]*(b/2) + * (qg[i] + qg[i+1]) - weight), + Mgust[i] >= (Mgust[i+1] + 0.5*deta[i]*(b/2) + * (Sgust[i] + Sgust[i+1])), + Mg[i] >= (Mg[i+1] + 0.5*deta[i]*(b/2) + * (Sg[i] + Sg[i+1])) + ]) + else: + constraints.extend([ + Sgust[i] >= (Sgust[i+1] + 0.5*deta[i]*(b/2) + * (qgust[i] + qgust[i+1])), + Sg[i] >= Sg[i+1] + 0.5*deta[i]*(b/2)*(qg[i] + qg[i+1]), + Mgust[i] >= (Mgust[i+1] + 0.5*deta[i]*(b/2) + * (Sgust[i] + Sgust[i+1])), + Mg[i] >= Mg[i+1] + 0.5*deta[i]*(b/2)*(Sg[i] + Sg[i+1]) + ]) self.submodels = [self.fs, self.aircraftPerf, self.slf, self.loading] @@ -643,21 +662,31 @@ def setup(self, aircraft, latitude=range(1, 21, 1), day=355): def test(): " test model for continuous integration " v = Aircraft(sp=False) - m = Mission(v, latitude=[15]) + m = Mission(v, latitude=[20]) m.cost = m[m.aircraft.Wtotal] m.solve() v = Aircraft(sp=True) - m = Mission(v, latitude=[15]) + m = Mission(v, latitude=[20]) m.cost = m[m.aircraft.Wtotal] m.localsolve() - v = Aircraft(Npod=5, sp=True) - m = Mission(v, latitude=[15]) + v = Aircraft(Npod=3, sp=True) + m = Mission(v, latitude=[20]) m.cost = m[m.aircraft.Wtotal] - m.localsolve() + f = relaxed_constants(M) + s = f.localsolve("mosek") + post_process(s) if __name__ == "__main__": SP = True - Vehicle = Aircraft(Npod=5, sp=SP) - M = Mission(Vehicle, latitude=[15]) + Vehicle = Aircraft(Npod=3, sp=SP) + M = Mission(Vehicle, latitude=[20]) M.cost = M[M.aircraft.Wtotal] - sol = (M.localsolve("mosek") if SP else M.solve("mosek")) + try: + sol = (M.localsolve("mosek") if SP else M.solve("mosek")) + except RuntimeWarning: + V2 = Aircraft(Npod=3, sp=SP) + M2 = Mission(V2, latitude=[20]) + M2.cost = M2[M2.aircraft.Wtotal] + feas = relaxed_constants(M2) + sol = feas.localsolve("mosek") + vks = post_process(sol)