How to use the gpkit.SignomialEquality function in gpkit

])
        
        #Horizontal Tail Constraints
        with SignomialsEnabled():
            constraints.extend([

                # Trim condition for each flight segment
                TCS([cruise['x_{AC}']/aircraft.wing['mac'] <= aircraft.wing['c_{m_{w}}']/cruise['C_{L}'] + \
                     cruise['x_{CG}']/aircraft.wing['mac'] + aircraft.HT['V_{ht}']*(cruise['C_{L_{ht}}']/cruise['C_{L}'])]),
                TCS([climb['x_{AC}']/aircraft.wing['mac'] <= aircraft.wing['c_{m_{w}}']/climb['C_{L}'] + \
                     climb['x_{CG}']/aircraft.wing['mac'] + aircraft.HT['V_{ht}']*(climb['C_{L_{ht}}']/climb['C_{L}'])]),


                aircraft.HT['L_{ht_{max}}'] == 0.5*aircraft['\\rho_0']*aircraft['V_{ne}']**2*aircraft.HT['S_{ht}']*aircraft.HT['C_{L_{ht,max}}'],
                #compute mrat, is a signomial equality
                SignomialEquality(aircraft.HT['m_{ratio}']*(1+2/aircraft.wing['AR']), 1 + 2/aircraft.HT['AR_{ht}']),

                #tail volume coefficient
                aircraft.HT['V_{ht}'] == aircraft.HT['S_{ht}']*aircraft.HT['x_{CG_{ht}}']/(aircraft.wing['S']*aircraft.wing['mac']),

                #enforce max tail location is the end of the fuselage
                aircraft.HT['x_{CG_{ht}}'] <= aircraft.fuse['l_{fuse}'],
                aircraft.HT['l_{ht}'] >= aircraft.HT['x_{CG_{ht}}'] - cruise['x_{CG}'],
                aircraft.HT['l_{ht}'] >= aircraft.HT['x_{CG_{ht}}'] - climb['x_{CG}'],

                #Stability constraint, is a signomial
                TCS([aircraft['SM_{min}'] + aircraft['\\Delta x_{CG}']/aircraft.wing['mac'] + aircraft.wing['c_{m_{w}}']/aircraft.wing['C_{L_{max}}'] <= aircraft.HT['V_{ht}']*aircraft.HT['m_{ratio}'] + aircraft.HT['V_{ht}']*aircraft.HT['C_{L_{ht,max}}']/aircraft.wing['C_{L_{max}}']]),

                # TCS([aircraft.wing['x_w'] >= cruise['x_{CG}'] + cruise['\\Delta x_w']]),
                # TCS([aircraft.wing['x_w'] >= climb['x_{CG}'] + climb['\\Delta x_w']]),

#Mhaero   >= rMh*Lhmax*(xtail-xwing), #[SP]
            #Mvaero   >= rMv*Lvmax*(xtail-xwing), #[SP]

            hfuse    == Rfuse, # may want to consider adding deltaRfuse later...
            xhbend   >= self.wingbox['x_{wing}'],

            
            # Horizontal bending model
            # Maximum axial stress is the sum of bending and pressurization stresses
            Ihshell <= ((pi+4*thetadb)*Rfuse**2)*Rfuse*tshell + 2/3*hdb**3*tdb, # [SP]
            #Ivshell <= (pi*Rfuse**2 + 8*wdb*Rfuse + (2*pi+4*thetadb)*wdb**2)*Rfuse*tshell, #[SP] #Ivshell approximation needs to be improved
            sigbend == rE*sigskin,
        
            # Horizontal bending material model
            # Calculating xbend, the location where additional bending material is required
            SignomialEquality(A0,A2*(xshell2-xhbend)**2 + A1*(xtail-xhbend)), #[SP] #[SPEquality] 
            A2      >=  self.ops['N_{land}']*(Wpay+Wshell+Wwindow+Winsul+Wfloor+Wseat)/(2*lshell*hfuse*sigMh), # Landing loads constant A2
            A1      >= (self.ops['N_{land}']*Wtail + rMh*self.htail['L_{h_{max}}'])/(hfuse*sigMh),                                # Aero loads constant A1
            A0      == (Ihshell/(rE*hfuse**2)),                                                # Shell inertia constant A0
            Ahbendf >= A2*(xshell2-self.wingbox['x_f'])**2 + A1*(xtail-self.wingbox['x_f']) - A0, #[SP]                           # Bending area forward of wingbox
            Ahbendb >= A2*(xshell2-self.wingbox['x_b'])**2 + A1*(xtail-self.wingbox['x_b']) - A0, #[SP]                           # Bending area behind wingbox

            Vhbendf >= A2/3*((xshell2-self.wingbox['x_f'])**3 - (xshell2-xhbend)**3) \
                            + A1/2*((xtail-self.wingbox['x_f'])**2 - (xtail - xhbend)**2) \
                            + A0*(xhbend-self.wingbox['x_f']), #[SP]

            Vhbendb >= A2/3*((xshell2-self.wingbox['x_b'])**3 - (xshell2-xhbend)**3) \
                            + A1/2*((xtail-self.wingbox['x_b'])**2 - (xtail - xhbend)**2) \
                            + A0*(xhbend-self.wingbox['x_b']), #[SP]
            Vhbendc >= .5*(Ahbendf + Ahbendb)*self.wingbox['c_0']*self.wingbox['\\bar_w'],
            Vhbend  >= Vhbendc + Vhbendf + Vhbendb,
            Whbend  >= g*rhobend*Vhbend,

#loose constraints on speed needed to prevent N from sliding out
                #to zero or infinity
                N1 >= .1,
                N1 <= 2,

                #note residuals 2 and 3 differ from TASOPT, by replacing mhc with mlc
                #in residual 4 I was able to remove the LPC/HPC mass flow equality
                #in residual 6 which allows for convergence
                #residual 2 HPT mass flow
                TCS([mhtD == (fp1)*mhc*(Pt25/Pt41)*(Tt41/Tt25)**.5]),
                
                #residual 3 LPT mass flow
                TCS([(fp1)*mlc*(Pt18/Pt45)*(Tt45/Tt18)**.5 == mltD]),
                
                #residual 4
                SignomialEquality(u7**2 +2*Cpfanex*T7, 2*Cpfanex*Tt7),
                rho7 == P7/(R*T7),
                TCS([mf*(Pt2/Pref)*(Tref/Tt2)**.5 == rho7*A7*u7]),
                
                #residual 5 core nozzle mass flow
                SignomialEquality(u5**2 +2*Cptex*T5, 2*Cptex*Tt5),
                rho5 == P5/(R*T5),
                
                #compute core mass flux
                mCore == rho5 * A5 * u5/(fp1),

                #compute fan mas flow
                mFan == rho7*A7*u7,
                
                #residual 6 LPC/HPC mass flow constraint
                TCS([mlc*(Pt18/Pref)*(Tref/Tt18)**.5 == mCore]),

#heat of combustion of jet fuel
        hf = Variable('h_f', 42.8, 'MJ/kg', 'Heat of Combustion of Jet Fuel')     #http://hypeRbook.com/facts/2003/EvelynGofman.shtml...prob need a better source

        #engine weight
        W_engine = Variable('W_{engine}', 'N', 'Weight of a Single Turbofan Engine')

        pilc = Variable('\pi_{lc}', '-', 'LPC Pressure Ratio')
        pihc = Variable('\pi_{hc}', '-', 'HPC Pressure Ratio')

        alpha = Variable('alpha', '-', 'By Pass Ratio')
        
        with SignomialsEnabled():
            constraints = [
                #making f+1 GP compatible --> needed for convergence
                SignomialEquality(fp1,f+1),

                #mass flow sizing
                mCore == Fd/(Fsp*a0*(alphap1)),  #B.194

                #component area sizing
                #fan area
                P2 == Pt2*(hold2)**(-3.512),
                T2 == Tt2 * hold2**-1,
                h2 == Cp1 * T2,
                rho2 == P2/(R * T2),  #B.196
                u2 == M2*(Cp1*R*T2/(781*units('J/kg/K')))**.5,  #B.197
                A2 == (alphap1)*mCore/(rho2*u2),     #B.198

                #HPC area
                P25 == Pt25*(hold25)**(-3.824857),
                T25 == Tt25 * hold25**-1,

sigth == dPover * Rfuse / tskin,

                # Floor loading
                lfloor >= lshell + 2 * Rfuse,
                Pfloor >= Nland * (Wpay + Wseat),
                Mfloor == 9. / 256. * Pfloor * wfloor,
                Afloor >= 2. * Mfloor / (sigfloor * hfloor) + 1.5 * Sfloor / taufloor,
                Vfloor == 2 * wfloor * Afloor,
                Wfloor >= rhofloor * g * Vfloor + 2 * wfloor * lfloor * Wppfloor,
                Sfloor == (5. / 16.) * Pfloor,
                hfloor <= 0.1 * Rfuse,

                # Tail cone sizing
                taucone == sigskin,
                Wcone >= rhocone * g * Vcone * (1 + fstring + fframe),
                SignomialEquality(xtail, lnose + lshell + .5 * lcone),  #[SP] #[SPEquality]
                
                # BENDING MODEL
                # Maximum axial stress is the sum of bending and pressurization
                # stresses
                Ihshell <= ((pi + 4 * thetadb) * Rfuse**2 + 8.*(1-thetadb**2/2) * (dRfuse/2.)*Rfuse + \
                            (2*pi + 4 * thetadb)*(dRfuse/2)**2) * Rfuse * tshell + \
                    2 / 3 * (hdb + dRfuse/2.)**3 * tdb,  # [SP]
                Ivshell <= (pi*Rfuse**2 + 8*wdb*Rfuse + (2*pi+4*thetadb)*wdb**2)*Rfuse*tshell, #[SP] #Ivshell
                # approximation needs to be improved
                sigbend == rE * sigskin,

                # Horizontal bending material model
                # Calculating xhbend, the location where additional bending
                # material is required
                xhbend >= xwing,
                SignomialEquality(A0h, A2h * (xshell2 - xhbend) ** 2 + A1h * (xtail - xhbend)), # [SP] #[SPEquality]

Vhbendf >= A2h / 3 * ((xshell2 - xf)**3 - (xshell2 - xhbend)**3) \
                + A1h / 2 * ((xtail - xf)**2 - (xtail - xhbend)**2) \
                + A0h * (xhbend - xf),  # [SP]

                Vhbendb >= A2h / 3 * ((xshell2 - xb)**3 - (xshell2 - xhbend)**3) \
                + A1h / 2 * ((xtail - xb)**2 - (xtail - xhbend)**2) \
                + A0h * (xhbend - xb),  # [SP]
                Vhbendc >= .5 * (Ahbendf + Ahbendb) * c0 * w,
                Vhbend >= Vhbendc + Vhbendf + Vhbendb,
                Whbend >= g * rhobend * Vhbend,

                # Wing variable substitutions
                SignomialEquality(xf, xwing + .5 * c0 * \
                                  w),  # [SP] [SPEquality]
                SignomialEquality(xb, xwing - .5 * c0 * \
                                  w),  # [SP] [SPEquality]

                sigMh <= sigbend - rE * dPover / 2 * Rfuse / tshell,
                SignomialEquality(xwing, lnose + 0.6 * lshell),  # TODO remove


                # Volume relations
                Vcyl == Askin * lshell,
                Vnose == Snose * tskin,
                Vbulk == Sbulk * tskin,
                Vdb == Adb * lshell,
                # TODO Revert to posynomial after debugging
                # [SP] #[SPEquality]
                SignomialEquality(
                    Vcabin, Afuse * (lshell + 0.67 * lnose + 0.67 * Rfuse)),

Pcabin = Variable('P_{cabin}','Pa','Cabin Air Pressure')
        W_buoy = Variable('W_{buoy}','lbf','Buoyancy Weight')
        Tcabin = Variable('T_{cabin}','K','Cabin Air Temperature')
        rhocabin = Variable('\\rho_{cabin}','kg/m^3','Cabin Air Density')

        constraints = []

        with SignomialsEnabled():
            constraints.extend([
            # Cabin Air properties
            rhocabin == Pcabin/(state['R']*Tcabin),
            Pcabin == 75000*units('Pa'),
            Tcabin == 297*units('K'),

            # Buoyancy weight #TODO relax the equality
            SignomialEquality(W_buoy,(rhocabin - state['\\rho'])*g*aircraft['V_{cabin}']),  #[SP] #[SPEquality]
            # W_buoy >= (rhocabin - state['\\rho'])*g*aircraft['V_{cabin}'], # [SP]

            # speed must be greater than stall speed
            state['V'] >= Vstall,

            # compute the drag
            D >= self.wingP['D_{wing}'] + self.fuseP['D_{fuse}'] + self.VTP['D_{vt}'] + self.HTP['D_{ht}'],
            LoD == W_avg/D,

            # Wing looading
            WLoad == .5 * self.wingP['C_{L}'] * self.aircraft['S'] * state.atm['\\rho'] * state['V']**2 / self.aircraft.wing['S'],

            # Geometric average of start and end weights of flight segment
            W_avg >= (W_start * W_end)**.5 + W_buoy, # Buoyancy weight included in Breguet Range

            # Maximum wing loading constraint

Pcabin = Variable('P_{cabin}','Pa','Cabin Air Pressure')
        W_buoy = Variable('W_{buoy}','lbf','Buoyancy Weight')
        Tcabin = Variable('T_{cabin}','K','Cabin Air Temperature')
        rhocabin = Variable('\\rho_{cabin}','kg/m^3','Cabin Air Density')

        constraints = []

        with SignomialsEnabled():
            constraints.extend([
            # Cabin Air properties
            rhocabin == Pcabin/(state['R']*Tcabin),
            Pcabin == 75000*units('Pa'),
            Tcabin == 297*units('K'),

            # Buoyancy weight #TODO relax the equality
            SignomialEquality(W_buoy,(rhocabin - state['\\rho'])*g*aircraft['V_{cabin}']),  #[SP] #[SPEquality]
            # W_buoy >= (rhocabin - state['\\rho'])*g*aircraft['V_{cabin}'], # [SP]

            # speed must be greater than stall speed
            state['V'] >= Vstall,

            # compute the drag
            D >= self.wingP['D_{wing}'] + self.fuseP['D_{fuse}'] + self.VTP['D_{vt}'] + self.HTP['D_{ht}'],
            LoD == W_avg/D,

            # Wing looading
            WLoad == .5 * self.wingP['C_{L}'] * self.aircraft['S'] * state.atm['\\rho'] * state['V']**2 / self.aircraft.wing['S'],

            # Geometric average of start and end weights of flight segment
            W_avg >= (W_start * W_end)**.5 + W_buoy, # Buoyancy weight included in Breguet Range

            # Maximum wing loading constraint