use relaxed constants to solve pod equations (#31)

* testing

* clean shear equations and solve with relaxed constants

* fix sens_chart and update solar for fixes in gplibrary

Author: Michael Burton

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,15 +78,15 @@ def plot_pods(df):
 
 def test():
     " for unit testing "
-    pods(N=[0])
+    pods(Nplot=100)
 
 if __name__ == "__main__":
     if len(sys.argv) > 1:
         path = sys.argv[1]
     else:
         path = ""
 
-    GENERATE = False
+    GENERATE = True
 
     if GENERATE:
         DF = pods()

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
+

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")

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)