How to use the gpkit.constraints.tight.TightConstraintSet function in gpkit

# Define vertical tail geometry

                Lvmax == 0.5*rho0*Vne**2*Svt*CLvmax,
                #compute the max force

                S == Svt*2,
                b == bvt*2,
                Lmax == 2*Lvmax,
                # Relate vertical tail geometry/weight to generic
                # wing used in structural model
     
                taper >= 0.27,
                # TODO: Constrain taper by tip Reynolds number
                # source: b737.org.uk

                TCS([xCGvt >= self.fuse['x_{CG}']+(dxlead+dxtrail)/2],
                         raiseerror=False),
                xCGvt <= self.fuse['l_{fuse}'],

                tau == tau,
                ])

        Model.__init__(self, None, constraints)

W_fcruise = Variable('W_{f_{cruise}}', 'N', 'Fuel Weight Burned in Cruise')
        W_total = Variable('W_{total}', 'N', 'Total Aircraft Weight')
        CruiseAlt = Variable('CruiseAlt', 'ft', 'Cruise Altitude [feet]')
        ReqRng = Variable('ReqRng', 'nautical_miles', 'Required Cruise Range')

        h = cls.state['h']
        hftClimb = cls.state['hft']
        dhft = cls.climbP['dhft']
        hftCruise = crs.state['hft']

        #make overall constraints
        constraints = []

        constraints.extend([
            #weight constraints
            TCS([ac['W_{e}'] + ac['W_{pay}'] + W_ftotal + ac['numeng'] * ac['W_{engine}'] + ac['W_{wing}'] <= W_total]),

            cls.climbP.aircraftP['W_{start}'][0] == W_total,
            cls.climbP.aircraftP['W_{end}'][-1] == crs.cruiseP.aircraftP['W_{start}'][0],

            # similar constraint 1
            TCS([cls.climbP.aircraftP['W_{start}'] >= cls.climbP.aircraftP['W_{end}'] + cls.climbP.aircraftP['W_{burn}']]),
            # similar constraint 2
            TCS([crs.cruiseP.aircraftP['W_{start}'] >= crs.cruiseP.aircraftP['W_{end}'] + crs.cruiseP.aircraftP['W_{burn}']]),

            cls.climbP.aircraftP['W_{start}'][1:] == cls.climbP.aircraftP['W_{end}'][:-1],
            crs.cruiseP.aircraftP['W_{start}'][1:] == crs.cruiseP.aircraftP['W_{end}'][:-1],

            TCS([ac['W_{e}'] + ac['W_{pay}'] + ac['numeng'] * ac['W_{engine}'] + ac['W_{wing}'] <= crs.cruiseP.aircraftP['W_{end}'][-1]]),

            TCS([W_ftotal >=  W_fclimb + W_fcruise]),
            TCS([W_fclimb >= sum(cls.climbP['W_{burn}'])]),

##                    TCS([ReqRngCruise   >= sum(RngCruise)]),
                    TCS([dhClimb2 + alt10k >= htoc])
                    ])
            for i in range(0, Nclimb1):
                constraints.extend([
                    SignomialEquality(hft[i], 1500*units('ft')+(i+1)*dhft[i])
                    ])

            for i in range(0, Nclimb2):
                constraints.extend([
                    SignomialEquality(hft[i+Nclimb1], 10000*units('ft')+(i+1)*dhft[i+Nclimb1])
                    ])

            for i in range(0, Nseg):
                constraints.extend([
                    TCS([W_start[i] >= W_end[i] + W_fuel[i]])
                    ])

        #constrain the segment weights in a loop
        for i in range(1, Nseg):
            constraints.extend([
                W_start[i] == W_end[i-1] 
                ])
        
        Model.__init__(self, W_ftotal, constraints, **kwargs)

constraints = []

        with SignomialsEnabled():
            #non vectorized constraints
            constraints.extend([
                LvtEO*lvt >= Te*y_eng + Dwm*y_eng,
                # Force moment balance for one engine out condition
                # TASOPT 2.0 p45

                TCS([dxlead + zmac*tanL + 0.25*cma >= lvt]), # [SP]
                # Tail moment arm

                LvtEO == 0.5*rho0*V1**2*Svt*CLvtEO,
                # Vertical tail force (y-direction) for engine out

                TCS([CLvtEO*(1 + clvtEO/(np.pi*e*Avt)) <= clvtEO]),
                #engine out CL computation

                Avt == bvt**2/Svt,

                Dwm >= 0.5*rho0*V1**2*self.engine['A_2']*CDwm,
                # Drag of a windmilling engine
                          
                Svt <= bvt*(croot + ctip)/2, # [SP]
                # Tail geometry relationship

                TCS([dxtrail >= croot + dxlead]),
                # Tail geometry constraint

                self.fuse['l_{fuse}'] >= dxtrail + self.fuse['x_{CG}'],
                # Fuselage length constrains the tail trailing edge

W_fcruise = Variable('W_{f_{cruise}}', 'N', 'Fuel Weight Burned in Cruise')
        W_total = Variable('W_{total}', 'N', 'Total Aircraft Weight')
        CruiseAlt = Variable('CruiseAlt', 'ft', 'Cruise Altitude [feet]')
        ReqRng = Variable('ReqRng', 'nautical_miles', 'Required Cruise Range')

        h = cls.state['h']
        hftClimb = cls.state['hft']
        dhft = cls.climbP['dhft']
        hftCruise = crs.state['hft']

        #make overall constraints
        constraints = []

        constraints.extend([
            #weight constraints
            TCS([ac['W_{e}'] + ac['W_{payload}'] + W_ftotal + ac['numeng'] * ac['W_{engine}'] + ac['W_{wing}'] + ac.VT['W_{struct}'] <= W_total]),

            cls.climbP.aircraftP['W_{start}'][0] == W_total,
            cls.climbP.aircraftP['W_{end}'][-1] == crs.cruiseP.aircraftP['W_{start}'][0],

            # similar constraint 1
            TCS([cls.climbP.aircraftP['W_{start}'] >= cls.climbP.aircraftP['W_{end}'] + cls.climbP.aircraftP['W_{burn}']]),
            # similar constraint 2
            TCS([crs.cruiseP.aircraftP['W_{start}'] >= crs.cruiseP.aircraftP['W_{end}'] + crs.cruiseP.aircraftP['W_{burn}']]),

            cls.climbP.aircraftP['W_{start}'][1:] == cls.climbP.aircraftP['W_{end}'][:-1],
            crs.cruiseP.aircraftP['W_{start}'][1:] == crs.cruiseP.aircraftP['W_{end}'][:-1],

            TCS([ac['W_{e}'] + ac['W_{payload}'] + ac['numeng'] * ac['W_{engine}'] + ac['W_{wing}'] + ac.VT['W_{struct}'] <= crs.cruiseP.aircraftP['W_{end}'][-1]]),

            TCS([W_ftotal >=  W_fclimb + W_fcruise]),
            TCS([W_fclimb >= sum(cls.climbP['W_{burn}'])]),

def __init__(self, **kwargs):
        T_tp = 216.65
        k = GRAVITATIONAL_ACCEL/(GAS_CONSTANT*T_tp)

        p11 = Variable('p_{11}', 22630, 'Pa', 'Pressure at 11 km')

        objective = 1/rho  # maximize density
        constraints = [h >= 11*units.km,
                       h <= 20*units.km,

                       # Temperature is constant in the tropopause
                       T == T_tp,

                       # Pressure-altitude relation, using taylor series exp
                       TCS([np.exp(k*11000)*p11/p >=
                            1 + te_exp_minus1(g/(R*T)*h, 15)], reltol=1E-4),

                       # Ideal gas law
                       rho == p/(R*T),
                       ]
        su = Sutherland()
        lc = su.link(constraints)

        Model.__init__(self, objective, lc, **kwargs)

if isinstance(surface, HorizontalTailNoStruct):
            AR = surface['AR_h']
            b = surface['b_{ht}']
            S = surface['S_h']
            p = surface['p_{ht}']
            q = surface['q_{ht}']
            tau = surface['\\tau_h']
            Lmax = surface['L_{{max}_h}']

        constraints = [
                       # Aspect ratio definition
                       AR == b**2/S,

                       # Defining taper dummy variables
                       TCS([p >= 1 + 2*taper]),
                       TCS([2*q >= 1 + p]),

                       # Upper bound on maximum thickness
                       tau <= 0.15,

                       # Root moment calculation (see Hoburg 2014)
                       # Depends on a given load the wing must support, Lmax
                       # Assumes lift per unit span proportional to local chord
                       Mr >= Lmax*AR*p/24,

                       # Root stiffness (see Hoburg 2014)
                       # Assumes rh = 0.75, so that rms box height = ~0.92*tmax
                       0.92*w*tau*tcap**2 + Icap <= 0.92**2/2*w*tau**2*tcap,

                       # Stress limit
                       # Assumes bending stress carried by caps (Icap >> Iweb)
                       8 >= Nlift*Mr*AR*q**2*tau/(S*Icap*sigmax),

w       = Variable('w', 0.5, '-', 'Wingbox-width-to-chord ratio')
        #xw     = Variable('x_w', 'm', 'Position of wing aerodynamic center')
        ymac    = Variable('y_{\\bar{c}_w}', 'm',
                           'Spanwise location of mean aerodynamic chord')


        objective = D

        with SignomialsEnabled():
            constraints = [
                           Lw == 0.5*rho*Vinf**2*Sw*CLw,

                           p >= 1 + 2*taper,
                           2*q >= 1 + p,
                           ymac == (b/3)*q/p,
                           TCS([(2./3)*(1+taper+taper**2)*croot/q <= cwma],
                               reltol=1E-2),
                           taper == ctip/croot,
                           TCS([Sw <= b*(croot + ctip)/2], reltol=1E-2), # [SP]

                           # DATCOM formula (Mach number makes it SP)
                           TCS([(AR/eta)**2*(1 + tanL**2 - M**2) + 8*pi*AR/CLaw
                                <= (2*pi*AR/CLaw)**2]),
                           CLw == CLaw*alpha,
                           alpha <= amax,

                           # Drag
                           D == 0.5*rho*Vinf**2*Sw*CDw,
                           CDw >= CDp + CLw**2/(pi*e*AR),
                           Re == rho*Vinf*cwma/mu,
                           1 >= (2.56*CLw**5.88/(Re**1.54*tau**3.32*CDp**2.62)
                              + 3.8e-9*tau**6.23/(CLw**0.92*Re**1.38*CDp**9.57)

#Breguet Range parameter constraints
                TCS([W_fuel[i]/W_end[i] >= te_exp_minus1(z_bre[i],3)]),
                TCS([z_bre[i] >= TSFC * t[i] * D[i]/ W_avg[i]]),

                #compute the Mach number
                M == V/a,

                #constraint on the segment time
                t[i] * V[i] == Rng[i],

                #total range is equal to the sum of the segment ranges
                TCS([ReqRng <= sum(Rng)]),

                #total fuel burn is equal to the sum of segment fuel burns
                TCS([W_fuelTotal <= sum(W_fuel)]),
                ]

        for j in range(1,n):
            constraints.extend([
                W_start[j] == W_end[j-1],
                ])

        #The objective is to minimze total fuel burn
        objective = W_fuelTotal  # Minimze total fuel burn

        #build the model
        Model.__init__(self, objective, constraints)

WLoad = Variable('W_{Load}', 'N/m^2', 'Wing Loading')
        t = Variable('tmin', 'min', 'Segment Flight Time in Minutes')
        thours = Variable('thr', 'hour', 'Segment Flight Time in Hours')

        constraints = []
        with SignomialsEnabled():
            constraints.extend([
                #speed must be greater than stall speed
                state['V'] >= Vstall,

                #Figure out how to delete
                Vstall == 120*units('kts'),
                WLoadmax == 6664 * units('N/m^2'),

                #compute the drag
                TCS([D >= self.wingP['D_{wing}'] + self.fuseP['D_{fuse}'] + self.htP['D_{ht}']]),

                #constraint CL and compute the wing loading
                W_avg == .5*self.wingP['C_{L}']*self.aircraft['S']*state.atm['\\rho']*state['V']**2,      
                WLoad == .5*self.wingP['C_{L}']*self.aircraft['S']*state.atm['\\rho']*state['V']**2/self.aircraft.wing['S'],

                #set average weight equal to the geometric avg of start and end weight
                W_avg == (W_start * W_end)**.5,

                #constrain the max wing loading
                WLoad <= WLoadmax,

                #compute fuel burn from TSFC
                W_burn == aircraft['numeng']*self.engineP['TSFC'] * thours * self.engineP['thrust'],
                   
                #time unit conversion
                t == thours,