How to use the gpkit.VectorVariable function in gpkit

W_end[:-1] >= W_end[1:] + W_fuel[1:], 
                            W_end[-1] >= W_zfw,
                            W_airframe >= f_airframe*MTOW,
                            W_zfw >= W_pay + W_avionics + W_airframe])
        
        #----------------------------------------------------
        # Steady level flight model

        CD = VectorVariable(Nseg, 'C_D', 0.05, '-', 'Drag coefficient')
    	CL = VectorVariable(Nseg, 'C_L', '-', 'Lift coefficient')
        V = VectorVariable(Nseg,'V', [15,15,15], 'm/s','cruise speed')
        rho = VectorVariable(Nseg, r'\rho', [1.2, 0.7, 1.2], 'kg/m^3', 'air density')
        S = Variable('S', 190, 'ft^2', 'wing area')
        eta_prop = VectorVariable(Nseg, r'\eta_{prop}', [0.6, 0.8, 0.6], '-',
                                  'propulsive efficiency')
        P_shaft = VectorVariable(Nseg, 'P_{shaft}', 'hp', 'Shaft power')

        constraints.extend([P_shaft >= V*(W_end+W_begin)/2*CD/CL/eta_prop,
                            0.5*rho*CL*S*V**2 == (W_end*W_begin)**0.5])

        #----------------------------------------------------
        # Breguet Range
        z_bre = VectorVariable(Nseg, 'z_{bre}', '-', 'breguet coefficient')
        BSFC = Variable('BSFC', 0.5, 'lbf/hr/hp', 'brake specific fuel consumption')
        t = VectorVariable(Nseg-1, 't', [5,0.2], 'days', 'time on station')
        R = Variable('R', 100, 'nautical_miles', 'range to station')
        g = Variable('g', 9.81, 'm/s^2', 'Gravitational acceleration')

        constraints.extend([z_bre[1:] >= V[1:]*t*BSFC*CD[1:]/CL[1:]/eta_prop[1:],
                            z_bre[0] >= R*BSFC*CD[0]/CL[0]/eta_prop[0],
                            W_fuel/W_end >= te_exp_minus1(z_bre, 3)])

# Propulsive efficiency variation with different flight segments,
        # will change depending on propeller characteristics

        #----------------------------------------------------
        # Engine Model (DF35)

        W_eng = Variable('W_{eng}', 'lbf', 'engine weight')
        W_engtot = Variable('W_{eng-tot}', 'lbf', 'Installed engine weight')
        W_engref = Variable('W_{eng-ref}', 4.4107, 'lbf',
                            'Reference engine weight')
        FuelOilFrac = Variable('FuelOilFrac', 0.98, '-', 'Fuel-oil fraction')
        P_shaftref = Variable('P_{shaft-ref}', 2.295, 'hp',
                              'reference shaft power')
        BSFC_min = Variable('BSFC_{min}', 0.32, 'kg/kW/hr', 'Minimum BSFC')
        BSFC = VectorVariable(NSeg, 'BSFC', 'lb/hr/hp',
                              'brake specific fuel consumption')
        RPM_max = Variable('RPM_{max}', 9000, 'rpm', 'Maximum RPM')
        RPM = VectorVariable(NSeg, 'RPM', 'rpm', 'Engine operating RPM')
        P_shaftmaxMSL = Variable('P_{shaft-maxMSL}', 'hp',
                                 'Max shaft power at MSL')
        P_shaftmax = VectorVariable(NSeg, 'P_{shaft-max}', 'hp',
                                    'Max shaft power at altitude')
        Lfactor = VectorVariable(NSeg, 'L_factor', '-',
                                 'Max shaft power loss factor')
        P_avn = Variable('P_{avn}', 40, 'watts', 'avionics power')
        P_pay = Variable('P_{pay}', 10, 'watts', 'payload power')
        P_shafttot = VectorVariable(NSeg, 'P_{shaft-tot}', 'hp',
                                    'total power, including avionics')
        eta_alternator = Variable(r'\eta_{alternator}', 0.8, '-',
                                  'alternator efficiency')
        P_shafttotmaxV = Variable('P_{shaft-tot-maxV}', 'hp',

R = Variable('R', 287, 'J/kg/K', 'Gas Constant for Air')

        #Engine
        TSFCcr2 = VectorVariable(Ncruise2, 'TSFC_{cr2}', '1/hr', 'Thrust Specific Fuel Consumption During Cruise2')
        thrustcr2 = VectorVariable(Ncruise2, 'thrust_{cr2}', 'N', 'Thrust During Cruise Segment #2')
        
        #defined here for linking purposes
        T0 = Variable('T_0', 'K', 'Free Stream Stagnation Temperature')
        P0 = Variable('P_0', 'kPa', 'Free Stream Static Pressure')

        #parameter to make breguent range eqn gp compatible
        z_brec2 = VectorVariable(Ncruise2, 'z_{brec2}', '-', 'Breguet Parameter')

        #range variables
        ReqRngCruise = Variable('ReqRngCruise', 'miles', 'Required Cruise Range')
        RngCruise2 = VectorVariable(Ncruise2, 'RngCruise2', 'miles', 'Segment Range During Cruise2')
        
        #Weights
        W_fuelCruise2 = VectorVariable(Ncruise2, 'W_{fuelCruise2}', 'N', 'Segment Fuel Weight')
        W_avgCruise2 = VectorVariable(Ncruise2, 'W_{avgCruise2}', 'N', 'Geometric Average of Segment Start and End Weight')
        W_avgClimb2 = VectorVariable(Nclimb2, 'W_{avgClimb2}', 'N', 'Geometric Average of Segment Start and End Weight')
        W_endCruise2 = VectorVariable(Ncruise2, 'W_{endCruise2}', 'N', 'Segment End Weight')

        #Fuselage area
        A_fuse = Variable('A_{fuse}', 'm^2', 'Estimated Fuselage Area')

        #non-GPkit variables
        #cruise 2 lsit
        izbre = map(int, np.linspace(0, Ncruise2 - 1, Ncruise2))
        
        constraints = []

# Avionics model
        W_avionics = Variable('W_{avionics}', 8, 'lbf', 'Avionics weight')
        Vol_avionics = Variable('Vol_{avionics}', 0.125, 'ft^3',
                                'Avionics volume')

        # end of first segment weight + first segment fuel weight must be
        # greater than MTOW.  Each end of segment weight must be greater than
        # the next end of segment weight + the next segment fuel weight. The
        # last end segment weight must be greater than the zero fuel weight
        constraints.extend([MTOW >= W_end[0] + W_fuel[0],
                            W_end[:-1] >= W_end[1:] + W_fuel[1:],
                            W_end[-1] >= W_zfw])

        #----------------------------------------------------
        # Steady level flight model
        CD = VectorVariable(NSeg, 'C_D', '-', 'Drag coefficient')
        CL = VectorVariable(NSeg, 'C_L', '-', 'Lift coefficient')
        V = VectorVariable(NSeg, 'V', 'm/s', 'cruise speed')
        S = Variable('S', 25.63, 'ft^2', 'wing area')
        eta_prop = VectorVariable(NSeg, '\\eta_{prop}', '-',
                                  'propulsive efficiency')
        eta_propCruise = Variable('\\eta_{prop-cruise}', 0.6, '-',
                                  'propulsive efficiency in cruise')
        eta_propClimb = Variable('\\eta_{prop-climb}', 0.5, '-',
                                 'propulsive efficiency in climb')
        eta_propLoiter = Variable('\\eta_{prop-loiter}', 0.7, '-',
                                  'propulsive efficiency in loiter')
        P_shaft = VectorVariable(NSeg, 'P_{shaft}', 'hp', 'Shaft power')
        T = VectorVariable(NSeg, 'T', 'lbf', 'Thrust')

        # Climb model
        # Currently climb rate is being optimized to reduce fuel consumption.

icruise2 = 0

izbre = map(int, np.linspace(0, Ncruise - 1, Ncruise))

#---------------------------------------------------------------------
#Variable definitions that apply to multiple classes
#physical constants
g = Variable('g', 9.81, 'm/s^2', 'Gravitational Acceleration')
gamma = Variable('\gamma', 1.4, '-', 'Air Specific Heat Ratio')
R = Variable('R', 287, 'J/kg/K', 'Gas Constant for Air')

#air properties
a = VectorVariable(Nseg, 'a', 'm/s', 'Speed of Sound')
rho = VectorVariable(Nseg, '\rho', 'kg/m^3', 'Air Density')
p = VectorVariable(Nseg, 'p', 'kPa', 'Pressure')
mu = VectorVariable(Nseg, '\mu', 'kg/m/s', 'Air Kinematic Viscosity')
T = VectorVariable(Nseg, 'T', 'K', 'Air Temperature')

#altitude
if Nclimb != 0:
    dhft = VectorVariable(Nclimb, 'dhft', 'feet', 'Change in Altitude Per Climb Segment [feet]') 
h = VectorVariable(Nseg, 'h', 'm', 'Altitude [meters]')
hft = VectorVariable(Nseg, 'hft', 'feet', 'Altitude [feet]')
dhClimb1 = Variable('dh_{climb1}', 8500, 'feet', 'Total Altitude Change Required in Climb 1')
dhClimb2 = Variable('dh_{climb2}', 'feet', 'Total Altitude Change Required in Climb 2')
dhTakeoff = Variable('dh_{takeoff}', 1500, 'feet', 'Total Altitude Change During Takeoff Profile')

#time
tmin = VectorVariable(Nseg, 'tmin', 'min', 'Flight Time in Minutes')
thours = VectorVariable(Nseg, 'thr', 'hour', 'Flight Time in Hours')

#Range

def setup(self):

        # define number of segments
        Nseg = 3

        constraints = []

        #----------------------------------------------------
        # Weight and fuel model

        MTOW = Variable('MTOW', 'lbf', 'max take off weight')
        W_end = VectorVariable(Nseg, 'W_{end}', 'lbf', 'segment-end weight')
        W_fuel = VectorVariable(Nseg, 'W_{fuel}', 'lbf',
                                'segment-fuel weight') 
        W_zfw = Variable('W_{zfw}', 'lbf', 'Zero fuel weight')
        W_pay = Variable('W_{pay}',10,'lbf', 'Payload weight')
        W_avionics = Variable('W_{avionics}', 2, 'lbf', 'Avionics weight')
        f_airframe = Variable('f_{airframe}', 0.3, '-', 'airframe weight fraction')
        W_airframe = Variable('W_{airframe}', 'lbf', 'airframe weight')
        W_begin = W_end.left # define beginning of segment weight
        W_begin[0] = MTOW 

        # end of first segment weight + first segment fuel weight must be greater 
        # than MTOW.  Each end of segment weight must be greater than the next end
        # of segment weight + the next segment fuel weight. The last end segment
        # weight must be greater than the zero fuel weight
        constraints.extend([MTOW >= W_end[0] + W_fuel[0], 
                            W_end[:-1] >= W_end[1:] + W_fuel[1:],

#excess power during climb
        excessPclimb2 = VectorVariable(Nclimb2, 'Excess Power Climb2', 'W', 'Excess Power During Climb')

        #climb rate and angle
        RCClimb2 = VectorVariable(Nclimb2, 'RCClimb2', 'feet/min', 'Rate of Climb/Decent')
        thetaClimb2 = VectorVariable(Nclimb2, '\\thetaClimb2', '-', 'Aircraft Climb Angle')

        #time
        tminClimb2 = VectorVariable(Nclimb2, 'tminClimb2', 'min', 'Flight Time in Minutes')
        thoursClimb2 = VectorVariable(Nclimb2, 'thrClimb2', 'hour', 'Flight Time in Hours')

        #range
        RngClimb2 = VectorVariable(Nclimb2, 'RngClimb2', 'miles', 'Segment Range During Climb2')

        #velocitites and mach numbers
        VClimb2 = VectorVariable(Nclimb2, 'VClimb2', 'knots', 'Aircraft Flight Speed')
        MClimb2 = VectorVariable(Nclimb2, 'MClimb2', '-', 'Aircraft Mach Number')

        #altitude
        dhftClimb2 = VectorVariable(Nclimb2, 'dhftClimb2', 'feet', 'Change in Altitude Per Climb Segment [feet]') 
        dhClimb2 = Variable('dh_{climb2}', 'feet', 'Total Altitude Change Required in Climb 2')

        W_fuelClimb2 = VectorVariable(Nclimb2, 'W_{fuelClimb2}', 'N', 'Segment Fuel Weight')
        W_avgClimb2 = VectorVariable(Nclimb2, 'W_{avgClimb2}', 'N', 'Geometric Average of Segment Start and End Weight')
        
        TSFCc2 = VectorVariable(Nclimb2, 'TSFC_{c2}', '1/hr', 'Thrust Specific Fuel Consumption During Climb2')
        thrustc2 = VectorVariable(Nclimb2, 'thrust_{c2}', 'N', 'Thrust During Climb Segment #2')

        thrustcr2 = VectorVariable(Ncruise2, 'thrust_{cr2}', 'N', 'Thrust During Cruise Segment #2')

        #Fuselage area
        A_fuse = Variable('A_{fuse}', 'm^2', 'Estimated Fuselage Area')

def __init__(self, N, **kwargs):

        #----------------------------------------------------
        # Fuel weight model

        W_start = Variable("W_{start}", "lbf",
                           "weight at beginning of flight segment")
        W_nplus1 = VectorVariable(N, "W_{N+1}", "lbf", "vector-end weight")
        W_fuel = VectorVariable(N, "W_{fuel}", "lbf",
                                "Segment-fuel weight")
        W_fuelfs = Variable("W_{fuel-fs}", "lbf",
                            "flight segment fuel weight")
        W_end = Variable("W_{end}", "lbf",
                         "weight at beginning of flight segment")
        W_n = VectorVariable(N, "W_{N}", "lbf", "vector-begin weight")

        # end of first segment weight + first segment fuel weight must be
        # greater  than MTOW.  Each end of segment weight must be greater
        # than the next end of segment weight + the next segment fuel weight.
        # The last end segment weight must be greater than the zero fuel
        # weight
        constraints = [W_start >= W_nplus1[0] + W_fuel[0],
                       W_fuelfs >= W_fuel.sum(),
                       W_nplus1[-1] >= W_end,
                       W_n[0] == W_start]

icruise2 = map(int, np.linspace(iclimb2[len(iclimb2)-1] + 1, Ncruise2 + iclimb2[len(iclimb2)-1], Ncruise2))
    ##    icruise2 = map(int, np.linspace(0, Ncruise2-1, Ncruise2))
    else:
        icruise2 = 0

    izbre = map(int, np.linspace(0, Ncruise - 1, Ncruise))

    #---------------------------------------------------------------------
    #Variable definitions that apply to multiple classes
    #physical constants
    g = Variable('g', 9.81, 'm/s^2', 'Gravitational Acceleration')
    gamma = Variable('\gamma', 1.4, '-', 'Air Specific Heat Ratio')
    R = Variable('R', 287, 'J/kg/K', 'Gas Constant for Air')

    #air properties
    a = VectorVariable(Nseg, 'a', 'm/s', 'Speed of Sound')
    rho = VectorVariable(Nseg, '\rho', 'kg/m^3', 'Air Density')
    p = VectorVariable(Nseg, 'p', 'kPa', 'Pressure')
    mu = VectorVariable(Nseg, '\mu', 'kg/m/s', 'Air Kinematic Viscosity')
    T = VectorVariable(Nseg, 'T', 'K', 'Air Temperature')

    #altitude
    if Nclimb != 0:
        dhft = VectorVariable(Nclimb, 'dhft', 'feet', 'Change in Altitude Per Climb Segment [feet]') 
    h = VectorVariable(Nseg, 'h', 'm', 'Altitude [meters]')
    hft = VectorVariable(Nseg, 'hft', 'feet', 'Altitude [feet]')
    dhClimb1 = Variable('dh_{climb1}', 8500, 'feet', 'Total Altitude Change Required in Climb 1')
    dhClimb2 = Variable('dh_{climb2}', 'feet', 'Total Altitude Change Required in Climb 2')
    dhTakeoff = Variable('dh_{takeoff}', 1500, 'feet', 'Total Altitude Change During Takeoff Profile')

    #time
    tmin = VectorVariable(Nseg, 'tmin', 'min', 'Flight Time in Minutes')

#defined here for linking purposes
        T0 = Variable('T_0', 'K', 'Free Stream Stagnation Temperature')
        P0 = Variable('P_0', 'kPa', 'Free Stream Static Pressure')

        #parameter to make breguent range eqn gp compatible
        z_bre = VectorVariable(Ncruise, 'z_{bre}', '-', 'Breguet Parameter')

        #range variables
        ReqRng = Variable('ReqRng', 'miles', 'Required Mission Range')
        RngCruise = VectorVariable(Ncruise, 'RngCruise', 'miles', 'Segment Range During Cruise')

        #Weight variables
        W_fuelCruise = VectorVariable(Ncruise, 'W_{fuelCruise}', 'N', 'Segment Fuel Weight')
        W_avgCruise = VectorVariable(Ncruise, 'W_{avgCruise}', 'N', 'Geometric Average of Segment Start and End Weight')
        W_avgClimb = VectorVariable(Nclimb, 'W_{avgClimb}', 'N', 'Geometric Average of Segment Start and End Weight')
        W_endCruise = VectorVariable(Ncruise, 'W_{endCruise}', 'N', 'Segment End Weight')

        #Fuselage area
        A_fuse = Variable('A_{fuse}', 'm^2', 'Estimated Fuselage Area')

        #non-GPkit variables
        #cruise 2 lsit
        icr = map(int, np.linspace(0, Ncruise - 1, Ncruise))

        constraints = []

        constraints.extend([
##            MCruise[icr] == 0.8,

            MCruise * aCruise == VCruise,