pgkirsch commented 9 years ago

I don't think SPs are returning sensitivities correctly. I have always been surprised by how few sensitivities get returned for some of my SPs. I'm aware that there is a threshold above which values must fall to get printed in results table, but I don't think I actually believe that so many of these fall below threshold. Consider two pieces of code (long again, sorry). I have an SP that only has one somewhat dumb signomial constraint in it at the moment. This signomial can easily be removed. The reason I haven't is a topic for another time. Anyway, if I solve the SP version, my results table has no sensitivities (surprising considering I have 30+ constants). But if I eliminate the signomial constraint thus making the almost exact same problem a GP and resolve, I get > 20 sensitivities. See below.

The original SP

from gpkit import Variable, Model, SignomialsEnabled, units
import numpy as np
import numpy.testing as npt

class Fuselage(Model):
    """
    Fuselage model
    """
    def setup(self):

        # Variables
#        Acone =    Variable('Acone', 'm^2', 'Cone enclosed area')
        Afloor =    Variable('Afloor', 'm^2', 'Floor area')
        Afuse =     Variable('Afuse', 'm^2', 'Cross sectional area')
        Askin =     Variable('Askin', 'm^2', 'Cross sectional area of skin')
        dRfuse =    Variable('dRfuse', 'm', 'Offset of lower fuselage cylinder')
        hfloor =    Variable('hfloor', 'm', 'Floor I-beam height') # Value
#        Ihshell =  Variable('Ihshell', 'm^4', 'Bending inertia')
#        Ivshell =  Variable('Ivshell', 'm^4', 'Bending inertia')
        lfloor =    Variable('lfloor', 'm', 'Floor length')
        lcone =     Variable('lcone', 'm', 'Cone length') # added this 
        lnose =     Variable('lnose', 'm', 'Nose length')
        lshell =    Variable('lshell', 'm', 'Shell length')
        Mfloor =    Variable('Mfloor', 'N*m', 'Maximum bending moment in floor beams')
        nrows =     Variable('n_{rows}', '-', 'Number of rows')
        Pfloor =    Variable('Pfloor', 'N', 'Distributed floor load')
        plamcone =  Variable('plamcone', '-', 'Cone taper ratio + 1')
        plamv =     Variable('plamv', '-', '1 + 2*Tail taper ratio')
        qlamv =     Variable('qlamv', '-', '1 + Tail taper ratio')
        Qv =        Variable('Qv', 'N*m', 'Torsion moment imparted by vertical tail')
#        Rcone =    Variable('Rcone', 'm', 'Cone radius')
        Rfuse =     Variable('Rfuse', 'm', 'Fuselage radius')
        Sbulk =     Variable('Sbulk', 'm^2', 'Bulkhead surface area')
        Sfloor =    Variable('Sfloor', 'N', 'Maximum shear in floor beams')
        sigth =     Variable('sigth', 'N/m^2', 'Hoop stress in skin')
        sigx =      Variable('sigx', 'N/m^2', 'Axial stress in skin')
        Snose =     Variable('Snose', 'm^2', 'Nose surface area')
        taucone =   Variable('taucone', 'N/m^2', 'Shear stress in cone')
        tcone =     Variable('tcone', 'm', 'Cone thickness')
        tshell =    Variable('tshell', 'm', 'Shell thickness')
        tskin =     Variable('tskin', 'm', 'Skin thickness')
        Vbulk =     Variable('Vbulk', 'm^3', 'Bulkhead skin volume')
        Vcone =     Variable('Vcone', 'm^3', 'Volume of cone skin')
        Vcyl =      Variable('Vcyl', 'm^3', 'Cylinder skin volume')
        Vfloor =    Variable('Vfloor', 'm^3', 'Floor volume')
        Vnose =     Variable('Vnose', 'm^3', 'Nose skin volume')
        Wapu =      Variable('Wapu', 'N', 'Weight of APU')
        Wcone =     Variable('Wcone', 'N', 'Weight of cone')
        wfloor =    Variable('wfloor', 'm', 'Floor width')
        Wfloor =    Variable('Wfloor', 'N', 'Floor weight')
        Wfuse =     Variable('Wfuse', 'N', 'Fuselage weight')
        Winsul =    Variable('Winsul', 'N', 'Weight of insulation material')
        Wpadd =     Variable('Wpadd', 'N', 'Misc weights (attendants, food, galley, toilets, doors etc.)')
        Wseat =     Variable('Wseat', 'N', 'Weight of seats')
        Wshell =    Variable('Wshell', 'N', 'Shell weight')
        Wskin =     Variable('Wskin', 'N', 'Skin weight')
        Wwindow =   Variable('Wwindow', 'N', 'Weight of windows')
#        xapu =     Variable('xapu', 120*0.3048, 'm', 'x-location of APU')
#        xconend =   Variable('xconend', 'm', 'x-location of cone end')
        xshell1 =   Variable('xshell1', 'm', 'Start of cylinder section') # value subbed as constraint
#        xshell2 =   Variable('xshell2', 102*0.3048, 'm', 'End of cyinder section')
#        xVbulk =   Variable('xVbulk', 'm^4', 'Volume moment of bulkhead')
#        xVcyl =    Variable('xVcyl', 'm^4', 'Volume moment of cylinder')
#        xVnose =   Variable('xVnose', 'm^4', 'Volume moment of nose')
#        xWapu =    Variable('xWapu', 'N*m', 'Moment of APU')
#        xWcone =   Variable('xWcone', 'N*m', 'Moment of cone')
#        xWfix =    Variable('xWfix', 'N*m', 'Moment of fixed weights')
#        xWfloor =  Variable('xWfloor', 'N*m', 'Moment of floor weight')
#        xWfuse =   Variable('xWfuse', 'N*m', 'Fuselage moment')
#        xWinsul =  Variable('xWinsul', 'N*m', 'Mass moment of insulation material')
#        xWpadd =   Variable('xWpadd', 'N*m', 'Moment of misc weights')
#        xWseat =   Variable('xWseat', 'N*m', 'Moment of seats')
#        xWshell =  Variable('xWshell', 'N*m', 'Mass moment of shell')
#        xWskin =   Variable('xWskin', 'N*m', 'Mass moment of skin')
#        xWwindow = Variable('xWwindow', 'N*m', 'Mass moment of windows')

        # Constants ([TAS] means sourced from TASOPT)
        bv =        Variable('bv', 7, 'm', 'Vertical tail span')
#        CLvmax =   Variable('CLvmax', 2.6, '-', 'Tail max lift coefficient') # [TAS]
        dp =        Variable('\\Deltap', 76000-24000, 'Pa', 'Pressure difference across fuselage skin')
        fapu =      Variable('fapu', 0.035, '-', 'APU weight as fraction of payload weight') # [TAS]
        ffadd =     Variable('ffadd', 0.20, '-', 'Skin proportional added weight of local reinforcmenets') # [TAS]
        fframe =    Variable('fframe', 0.25, '-', 'Fractional weight of frame')
        fpadd =     Variable('fpadd', 0.35, '-', 'All other misc weight as fraction of payload weight') # [TAS]
        fseat =     Variable('fseat', 0.10, '-', 'Seat weight as fraction of payload weight')
        fstring =   Variable('fstring', 0.35, '-', 'Fraction of stringers?') # [TAS]
        g =         Variable('g', 9.81, 'm/s^2', 'Gravitational acceleration')
        Lvmax =     Variable('Lvmax', 35000, 'N', 'Tail max lift')
        Nland =     Variable('Nland', 6.0,  '-', 'Emergency landing load factor') # [TAS]
        nseats =    Variable('n_{seats}', 186, '-', 'Number of seats') # Must be multiple of SPR
        pitch =     Variable('pitch', 31 , 'in', 'Seat pitch') # Boeing doc
#        qNE =      Variable('qNE', 'Pa', 'Never exceed dynamic pressure')
        rE =        Variable('rE', 1.0, '-', 'Ratio of stringer/skin moduli') # [TAS]
        rhobend =   Variable('rhobend', 2700, 'kg/m^3', 'Stringer density') # [TAS]
        rhocone =   Variable('rhocone', 2700, 'kg/m^3', 'Density of cone material') # [TAS]
        rhofloor =  Variable('rhofloor', 2700, 'kg/m^3', 'Density of floor material') # [TAS]
        rhoskin =   Variable('rhoskin', 2700, 'kg/m^3', 'Skin density') # [TAS]
        sigfloor =  Variable('sigfloor', 30000/0.000145, 'N/m^2', 'Max allowable cap stress') # [TAS]
        sigskin =   Variable('sigskin', 15000/0.000145, 'N/m^2', 'Max allowable skin stress') # [TAS]
        SPR =       Variable('SPR', 6, '-', 'Number of seats per row')
#        Sv =       Variable('Sv', 20, 'm^2', 'Tail reference area')
        taufloor =  Variable('taufloor', 30000/0.000145, 'N/m^2', 'Max allowable shear web stress') # [TAS]
        waisle =    Variable('w_{aisle}', 0.51, 'm', 'Aisle width') # Boeing doc
        Wfix =      Variable('Wfix', 3000*4.45, 'N', 'Fixed weights (pilots, cockpit seats, instrumentation, navcom)')
        Wpay =      Variable('Wpay', 177000, 'N', 'Payload weight') # b737.org.uk
        Wppfloor =  Variable('Wppfloor', 60, 'N/m^2', 'Floor weight/area density') # [TAS]
        Wppinsul =  Variable('Wppinsul', 22, 'N/m^2', 'Weight/Area density of insulation material') # [TAS]
        Wpwindow =  Variable('Wpwindow', 145.*3, 'N/m', 'Weight/Length density of windows') # [TAS]
        wseat =     Variable('w_{seat}', 0.5, 'm', 'Seat width') # Boeing doc
#        xfix =     Variable('xfix', 7*0.3048, 'm', 'x-location of fixed weights')

#        Vcabin
#        rhocabin
#        rho0
#        R
#        Tcabin
#        pcabin
#        Wbuoy
#        Weff

        with SignomialsEnabled():
            objective = Wfuse

            constraints = [
                            # Geometry relations
                            lnose == xshell1,

                            # Cross section relations
                            Rfuse >= 77*0.0254*units.m, # TODO: this needs to be driven by something e.g. number of passengers
                            dRfuse >= 15*0.0254*units.m, # TODO: this needs to be driven by something e.g. cargo capacity
                            Askin >= 2*np.pi*Rfuse*tskin + 2*dRfuse*tskin,
                            Afuse >= np.pi*Rfuse**2,# + 2*Rfuse*dRfuse, # TODO Fix this
                            tshell >= tskin*(1 + rE*fstring*rhoskin/rhobend),
    #                            Ihshell >= (np.pi*Rfuse**2 + 2*np.pi*(dRfuse/2)**2)*Rfuse*tshell,
    #                            Ivshell >= np.pi*Rfuse**2 * Rfuse*tshell,

                            # Pressure shell loads
                            sigx == dp/2*Rfuse/tshell,
                            sigth == dp*Rfuse/tskin,
                            tskin == dp*Rfuse/sigskin,
                            Snose**(8./5) >= (2*np.pi*Rfuse**2)**(8./5) * (1./3 + (2./3)*(lnose/Rfuse)**(8./5)),
                            Sbulk == 2*np.pi*Rfuse**2,
                            Vcyl == Askin*lshell,
                            Vnose == Snose*tskin,
                            Vbulk == Sbulk*tskin,
    #                        xVcyl >= 0.5*(xshell1+xshell2)*Vcyl,
    #                        xVnose >= 0.5*(xshell1)*Vnose,
    #                        xVbulk >= (xshell2 + 0.5*dRfuse)*Vbulk,
                            Wskin >= rhoskin*g*(Vcyl + Vnose + Vbulk),
    #                        xWskin >= rhoskin*g*(xVcyl + xVnose + xVbulk),
                            Wshell >= Wskin*(1 + fstring + fframe + ffadd),
    #                        xWshell >= xWskin*(1 + fstring + fframe + ffadd),

                            # Cabin volume and buoyancy weight
    #                            rhocabin == (1/R*Tcabin)*pcabin,
    #                            Wbuoy >= (rhocabin - rho0)*g*Vcabin, # SP
    #                            Weff >= W + Wbuoy,

                           # Windows and insulation
                           Wwindow == Wpwindow * lshell,
    #                       xWwindow >= 0.5* (xshell1+xshell2)*Wwindow,
                           Winsul >= Wppinsul*(1.1*np.pi*Rfuse*lshell + 0.55*(Snose + Sbulk)),
    #                       xWinsul >= 0.5*(xshell1 + xshell2)*Winsul,

                           # Payload-proportional weights
                           Wapu == Wpay*fapu,
    #                       xWapu == xapu*Wapu,
                           Wseat == Wpay*fseat,
    #                       xWseat >= 0.5*(xshell1 + xshell2)*Wseat,
                           Wpadd == Wpay*fpadd,
    #                       xWpadd >= 0.5*(xshell1 + xshell2)*Wpadd,

                            # Fixed weight
    #                       xWfix == xfix*Wfix,

                            # Floor
                            Pfloor >= Nland*(Wpay + Wseat),
                            Sfloor == 0.5*Pfloor, # without support
                            Mfloor == 0.25*Pfloor*wfloor,
#                            lshell >= xshell2 - xshell1, # NOTE: [SP]
#                            xshell2 >= xshell1 + lshell, # NOTE: requires downward pressure on xshell2
                            lnose >= 5.18*units.m, # TODO: this needs to be driven by something e.g cockpit length
                            Afloor >= 2*Mfloor/(sigfloor*hfloor) + 1.5*Sfloor/taufloor,
                            hfloor <= 5*0.0254*units.m, # TODO: Get better upper bound on hfloor
                            Vfloor >= 2*Rfuse*Afloor,
                            lfloor >= lshell + 2*Rfuse,
                            Wfloor >= rhofloor*g*Vfloor + 2*wfloor*lfloor*Wppfloor,
   #                        xWfloor >= 0.5*(xshell1 + xshell2)*Wfloor,

                            # Tail cone
                            Qv == Lvmax*bv*plamv/(3*qlamv), # uses same p, q trick as GP for A/C design
                            2*qlamv <= 1 + plamv, # NOTE: [SP] (could be eliminated if necessary)
                            plamv >= 1.6,
                            taucone == sigskin,
                            Qv == 2*Afuse*taucone*tcone, # modified to match TASOPT code
    #                           Acone >= Afuse*(1 + plamcone*qlamcone*()/())
    #                           Rcone >= ....
    #                           tcone
                            Vcone*Rfuse*plamcone >= 2*Qv/taucone * lcone,
                            lcone >= 15*0.3048*units.m, # TODO: this needs to be driven by something (e.g. APU size)
    #                       xshell2 + lcone <= xconend,
                            plamcone <= 1.3,
                            Wcone >= rhocone*g*Vcone*(1 + fstring + fframe + ffadd),
    #                        xWcone >= 0.5*(xshell2 + xconend) * Wcone,

                            # Passenger requirements
                            nseats == nrows*SPR,
                            lshell >= nrows*pitch,
                            wfloor >= SPR*wseat + waisle,

                            # Total fuselage weight
                            Wfuse >= Wfix + Wapu + Wpadd + Wseat + Wshell + Wwindow + Winsul
                                     + Wcone + Wfloor,
    #                        xWfuse >= xWfix + xWapu + xWpadd + xWseat + xWshell + xWcone
    #                                 + xWwindow + xWinsul + xWfloor
                          ]
            return objective, constraints

    def test(self):
        sol = self.localsolve()

#        npt.assert_almost_equal(sol('lshell'), sol('xshell2') - sol('xshell1'), decimal=2)
        npt.assert_almost_equal(2*sol('qlamv'), 1 + sol('plamv'), decimal=5)

if __name__ == "__main__":
    Fuselage().test()

################
$ ipython fuselage.py

....

Sensitivities
-------------