def solve_jiang_guan_lp(return_samples, targets, asset_partition):
    n = return_samples.shape[1]
    prob = pulp.LpProblem('Jiang Guan DCCP Sample Approximation', pulp.LpMinimize)
    x = [pulp.LpVariable('Asset {0:d} Weight'.format(i), 0.0, 1.0) for i in range(n)]
    mu = calc_first_moment(return_samples)
    prob += pulp.lpDot(mu, x), 'Sample Mean Return'

    i = 0
    for sample in return_samples:
        prob += (pulp.lpDot(sample, x) >= targets[0],
                 'Global return target for sample {0:d}'.format(i))
        i += 1

    i = 0
    j = 1
    for assets in asset_partition:
        for sample in return_samples:
            prob += (pulp.lpSum([sample[k] * x[k] for k in range(n)]) >= targets[j],
                     'Return target for segment {0:d} for sample {1:d}'.format(j, i))
            i += 1
        j += 1

    prob += (pulp.lpSum(x) == 1.0, 'Fully invested portfolio requirement')

    prob.writeLP('JiangGuanDccp.lp')
    prob.solve()
    status = pulp.LpStatus[prob.status]
    print 'Status: {0:s}'.format(status)
    return np.array([v.varValue for v in prob.variables()]), status

def solve(g):
    el = g.get_edge_list()
    nl = g.get_node_list()
    p = LpProblem('min_cost', LpMinimize)
    capacity = {}
    cost = {}
    demand = {}
    x = {}
    for e in el:
        capacity[e] = g.get_edge_attr(e[0], e[1], 'capacity')
        cost[e] = g.get_edge_attr(e[0], e[1], 'cost')
    for i in nl:
        demand[i] = g.get_node_attr(i, 'demand')
    for e in el:
        x[e] = LpVariable("x"+str(e), 0, capacity[e])
    # add obj
    objective = lpSum (cost[e]*x[e] for e in el)
    p += objective
    # add constraints
    for i in nl:
        out_neig = g.get_out_neighbors(i)
        in_neig = g.get_in_neighbors(i)
        p += lpSum(x[(i,j)] for j in out_neig) -\
             lpSum(x[(j,i)] for j in in_neig)==demand[i]
    p.solve()
    return x, value(objective)

    def _GetTotalEnergyProblem(self, min_driving_force=0, objective=pulp.LpMinimize):
        
        # Define and apply the constraints on the concentrations
        ln_conc_lb, ln_conc_ub = self._MakeLnConcentratonBounds()

        # Create the driving force variable and add the relevant constraints
        A, b, _c = self._MakeDrivingForceConstraints(ln_conc_lb, ln_conc_ub)
       
        lp = pulp.LpProblem("OBD", objective)
        
        # ln-concentration variables
        l = pulp.LpVariable.dicts("l", ["%d" % i for i in xrange(self.Nc)])
        x = [l["%d" % i] for i in xrange(self.Nc)] + [min_driving_force]
        
        total_g = pulp.LpVariable("g_tot")
        
        for j in xrange(A.shape[0]):
            row = [A[j, i] * x[i] for i in xrange(A.shape[1])]
            lp += (pulp.lpSum(row) <= b[j, 0]), "energy_%02d" % j
        
        total_g0 = float(self.dG0_r_prime * self.fluxes.T)
        total_reaction = self.S * self.fluxes.T
        row = [total_reaction[i, 0] * x[i] for i in xrange(self.Nc)]
        lp += (total_g == total_g0 + pulp.lpSum(row)), "Total G"

        lp.setObjective(total_g)
        
        #lp.writeLP("../res/total_g.lp")
        
        return lp, total_g

    def K_dominnace_check_2(self, u_d, v_d, _inequalities):
        """

        :param u_d: a d-dimensional vector(list) like [ 8.53149891  3.36436796]
        :param v_d: tha same list like u_d
        :param _inequalities: list of constraints on d-dimensional Lambda Polytope like
         [[0, 1, 0], [1, -1, 0], [0, 0, 1], [1, 0, -1], [0.0, 1.4770889, -3.1250839]]
        :return: True if u is Kdominance to v regarding given _inequalities otherwise False
        """
        _d = len(u_d)

        prob = LpProblem("Kdominance", LpMinimize)
        lambda_variables = LpVariable.dicts("l", range(_d), 0)

        for inequ in _inequalities:
            prob += lpSum([inequ[j + 1] * lambda_variables[j] for j in range(0, _d)]) + inequ[0] >= 0

        prob += lpSum([lambda_variables[i] * (u_d[i]-v_d[i]) for i in range(_d)])

        #prob.writeLP("show-Ldominance.lp")

        status = prob.solve()
        LpStatus[status]

        result = value(prob.objective)
        if result < 0:
            return False

        return True

 def _MakeMDFProblem(self):
     """Create a CVXOPT problem for finding the Maximal Thermodynamic
     Driving Force (MDF).
    
     Does not set the objective function... leaves that to the caller.
    
     Returns:
         the linear problem object, and the three types of variables as arrays
     """
     A, b, c, y, l = self._GetPrimalVariablesAndConstants()
     B = pulp.LpVariable("mdf")
     x = y + l + [B]
     lp = pulp.LpProblem("MDF_PRIMAL", pulp.LpMaximize)
     
     cnstr_names = ["driving_force_%02d" % j for j in xrange(self.Nr_active)] + \
                   ["covariance_var_ub_%02d" % j for j in xrange(self.Nr)] + \
                   ["covariance_var_lb_%02d" % j for j in xrange(self.Nr)] + \
                   ["log_conc_ub_%02d" % j for j in xrange(self.Nc)] + \
                   ["log_conc_lb_%02d" % j for j in xrange(self.Nc)]
       
     for j in xrange(A.shape[0]):
         row = [A[j, i] * x[i] for i in xrange(A.shape[1])]
         lp += (pulp.lpSum(row) <= b[j, 0]), cnstr_names[j]
     
     objective = pulp.lpSum([c[i] * x[i] for i in xrange(A.shape[1])])
     lp.setObjective(objective)
     
     lp.writeLP("res/mdf_primal.lp")
     
     return lp, objective, y, l, B

def good_annotation_locations(item, annotate_first=True):
    """Find the minimum number of annotations necessary to extract all the fields

    Since annotations can be reviewed and modified later by the user we want to keep
    just the minimum number of them.

    Parameters
    ----------
    item : Item
    annotate_first : book
        If true always annotate the first instance of the item in the page

    Returns
    -------
    List[ItemLocation]
    """
    #    x[i] = 1 iff i-th item is representative
    # A[i, j] = 1 iff i-th item contains the j-th field
    #
    # Solve:
    #           min np.sum(x)
    # Subject to:
    #           np.all(np.dot(A.T, x) >= np.repeat(1, len(fields)))
    index_locations = {location: i for i, location in enumerate(item.locations)}
    P = pulp.LpProblem("good_annotation_locations", pulp.LpMinimize)
    X = [pulp.LpVariable("x{0}".format(i), cat="Binary") for i in range(len(index_locations))]
    if annotate_first:
        P += X[0] == 1
    P += pulp.lpSum(X)
    for field in item.fields:
        P += pulp.lpSum([X[index_locations[location.item]] for location in field.locations]) >= 1
    P.solve()
    return [i for (i, x) in enumerate(X) if x.value() == 1]

    def _MakeMDFProblemDual(self):
        """Create a CVXOPT problem for finding the Maximal Thermodynamic
        Driving Force (MDF).
       
        Does not set the objective function... leaves that to the caller.
       
        Returns:
            the linear problem object, and the four types of variables as arrays
        """
        A, b, c, w, g, z, u = self._GetDualVariablesAndConstants()
        x = w + g + z + u
        lp = pulp.LpProblem("MDF_DUAL", pulp.LpMinimize)

        cnstr_names = ["y_%02d" % j for j in xrange(self.Nr)] + \
                      ["l_%02d" % j for j in xrange(self.Nc)] + \
                      ["MDF"]
        
        for i in xrange(A.shape[1]):
            row = [A[j, i] * x[j] for j in xrange(A.shape[0])]
            lp += (pulp.lpSum(row) == c[i, 0]), cnstr_names[i]

        objective = pulp.lpSum([b[i] * x[i] for i in xrange(A.shape[0])])
        lp.setObjective(objective)
        
        lp.writeLP("res/mdf_dual.lp")
        
        return lp, objective, w, g, z, u

def min_one_norm(B,initial_seed,seed):

    weight_initial = 1 / float(len(initial_seed))
    weight_later_added = weight_initial / float(0.5)
    difference = len(seed) - len(initial_seed)
    [r,c] = B.shape
    prob = pulp.LpProblem("Minimum one norm", pulp.LpMinimize)
    indices_y = range(0,r)
    y = pulp.LpVariable.dicts("y_s", indices_y, 0)
    indices_x = range(0,c)
    x = pulp.LpVariable.dicts("x_s", indices_x)

    f = dict(zip(indices_y, [1.0]*r))

    prob += pulp.lpSum(f[i] * y[i] for i in indices_y) # objective function
    
    prob += pulp.lpSum(y[s] for s in initial_seed) >= 1

    prob += pulp.lpSum(y[r] for r in seed) >= 1 + weight_later_added * difference

    for j in range(r):
        temp = dict(zip(indices_x, list(B[j,:])))
        prob += pulp.lpSum(y[j] + (temp[k] * x[k] for k in indices_x)) == 0

    prob.solve()

    print "Status:", pulp.LpStatus[prob.status]
    result = []
    for var in indices_y:
        result.append(y[var].value())
   
    return result 

def opt(C, X):
    orderNum = len(X)
    routeNum = len(C)
    routeIdx = range(routeNum)
    orderIdx = range(orderNum)
    # print routeIdx,orderIdx
    eps = 1.0 / 10 ** 7
    print eps
    var_choice = lp.LpVariable.dicts('route', routeIdx, cat='Binary')
    # var_choice=lp.LpVariable.dicts('route',routeIdx,lowBound=0)#尝试松弛掉01变量
    exceed_labor = lp.LpVariable('Number of routes exceed 1000', 0)
    prob = lp.LpProblem("lastMile", lp.LpMinimize)

    prob += exceed_labor * 100000 + lp.lpSum(var_choice[i] * C[i] for i in routeIdx)

    prob += lp.lpSum(var_choice[i] for i in routeIdx) <= 1000 + exceed_labor + eps
    for i in orderIdx:
        prob += lp.lpSum(var_choice[j] for j in X[i]) >= (1 - eps)

    prob.solve(lp.CPLEX(msg=0))
    print "\n\nstatus:", lp.LpStatus[prob.status]
    if lp.LpStatus[prob.status] != 'Infeasible':
        obj = lp.value(prob.objective)
        print "\n\nobjective:", obj
        sol_list = [var_choice[i].varValue for i in routeIdx]
        print "\n\nroutes:", (sum(sol_list))
        # print "\n\noriginal problem:\n",prob
        return obj, sol_list, lp.LpStatus[prob.status]
    else:
        return None, None, lp.LpStatus[prob.status]

    def fit(self, x, y):
        classifiers = []
        classifier_weights = []
        clf = self.base_estimator


        # get predictions of one classifer
        clf.fit(x, y)
        y_predict = clf.predict(x)
        u = (y_predict == y).astype(int)
        u[u == 0] = -1
        for index, value in enumerate(u):
            if value == -1:
                print index


        # solving linear programming
        d = pp.LpVariable.dicts("d", range(len(y)), 0, 1)
        prob = pp.LpProblem("LPboost", pp.LpMinimize)
        prob += pp.lpSum(d) == 1  # constraint for sum of weights to be 1
        # objective function
        objective_vector = []
        for index in range(len(y)):
            objective_vector.append(d[index] * u[index])
        prob += pp.lpSum(objective_vector)


        print pp.LpStatus[prob.solve()]
        for v in prob.variables():
            if v.varValue > 0:
                print v.name + "=" + str(v.varValue)

        def setup():
            # ... declare variables
            x = ilp.LpVariable.dicts('x', nodes_vars.keys(), 0, 1, ilp.LpBinary)
            y = ilp.LpVariable.dicts('y',
                                     [(i, j) for i, j in edges] + [(j, i) for i, j in edges],
                                     0, 1, ilp.LpBinary)
            limits = defaultdict(int)
            for i, j in edges:
                limits[i] += 1
                limits[j] += 1

            # ... define the problem
            prob = ilp.LpProblem("Factorizer", ilp.LpMinimize)

            # ... define the constraints
            for i in nodes_vars:
                prob += ilp.lpSum(y[(i, j)] for j in nodes_vars if (i, j) in y) <= limits[i]*x[i]

            for i, j in edges:
                prob += y[(i, j)] + y[(j, i)] == 1

            # ... define the objective function (min number of factorizations)
            prob += ilp.lpSum(x[i] for i in nodes_vars)

            return x, prob

    def K_dominance_check(self, _V_best_d, Q_d):
        """
        :param _V_best_d: a list of d-dimension
        :param Q_d: a list of d-dimension
        :return: True if _V_best_d is prefered to Q_d regarding self.Lambda_inequalities and using Kdominance
         other wise it returns False
        """
        _d = len(_V_best_d)

        prob = LpProblem("Ldominance", LpMinimize)
        lambda_variables = LpVariable.dicts("l", range(_d), 0)

        for inequ in self.Lambda_ineqalities:
            prob += lpSum([inequ[j + 1] * lambda_variables[j] for j in range(0, _d)]) + inequ[0] >= 0

        prob += lpSum([lambda_variables[i] * (_V_best_d[i]-Q_d[i]) for i in range(_d)])

        #prob.writeLP("show-Ldominance.lp")

        status = prob.solve()
        LpStatus[status]
        result = value(prob.objective)

        if result < 0:
            return False

        return True

    def _create_lp(self):
        ''' See base class.
        '''
        self.model._create_lp()
        self.lp_model_max_slack = self.model.lp_model.deepcopy()

        input_slack_vars = pulp.LpVariable.dicts(
            'input_slack', self.strongly_disposal_input_categories,
            0, None, pulp.LpContinuous)
        output_slack_vars = pulp.LpVariable.dicts(
            'output_slack', self.strongly_disposal_output_categories,
            0, None, pulp.LpContinuous)

        # change objective function
        self.lp_model_max_slack.sense = pulp.LpMaximize
        self.lp_model_max_slack.objective = (
            pulp.lpSum(list(input_slack_vars.values())) +
            pulp.lpSum(list(output_slack_vars.values())))

        # change constraints
        for input_category in self.strongly_disposal_input_categories:
            name = self.model._constraints[input_category]
            self.lp_model_max_slack.constraints[name].addterm(
                input_slack_vars[input_category], 1)
            self.lp_model_max_slack.constraints[name].sense = pulp.LpConstraintEQ

        for output_category in self.strongly_disposal_output_categories:
            name = self.model._constraints[output_category]
            self.lp_model_max_slack.constraints[name].addterm(
                output_slack_vars[output_category], -1)
            self.lp_model_max_slack.constraints[name].sense = pulp.LpConstraintEQ

	def get_linear_program_solution(self, c, b, A, x):
		prob = LpProblem("myProblem", LpMinimize)
		prob += lpSum([xp*cp for xp, cp in zip(x, c)]), "Total Cost of Ingredients per can" 

		for row, cell in zip(A, b):
			prob += lpSum(ap*xp for ap, xp in zip(row,  x)) <= cell

		solved = prob.solve()
		return prob

    def def_problem(self, queue_a):
        """define lp problem"""
        n = self.n
        tn = self.tn
        tl = self.tl
        b = self.b

        iIdx = range(n)
        tIdx = range(tn)
        ti = range(n)
        si = range(n)
        Ti = range(n)
        Si = range(n)

        for i in xrange(n):
            ti[i] = queue_a[i].at
            si[i] = queue_a[i].tsize
            Ti[i] = queue_a[i].ft
            Si[i] = queue_a[i].size
            iIdx[i] = str(i)

        for t in xrange(tn):
            tIdx[t] = str(t)

        """-----------------------------------
        PuLP variable definition
        varb--bi(t)
        prob--objective and constraints

        objective:
            max:[0~n-1][0~ti-1])sigma bi(t)

        constraints:
            1.any t,[0~n-1] sigma bi(t)   <= B
            2.any i,[0-Ti]  sigma bi(t)*tl>= si
            3.any i,[0~T-1] sigma bi(t)*tl<= Si
        ------------------------------------"""
        print "\ndefine lp variables"
        self.varb = pulp.LpVariable.dicts('b',(iIdx,tIdx),0,b,cat='Integer')

        print "define lp problem"
        self.prob = pulp.LpProblem('Prefetching Schedule',pulp.LpMaximize)

        print "define lp objective"
        self.prob += pulp.lpSum([tl*self.varb[i][t] for i in iIdx for t in tIdx if int(t)<ti[int(i)]])

        print "define constraints on B"
        for t in tIdx:
            self.prob += pulp.lpSum([self.varb[i][t] for i in iIdx]) <= b

        print "define constraints on si"
        for i in iIdx:
            self.prob += pulp.lpSum([tl*self.varb[i][t] for t in tIdx if int(t)<=Ti[int(i)]]) >= si[int(i)]

        print "define constraints on Si"
        for i in iIdx:
            self.prob += pulp.lpSum([tl*self.varb[i][t] for t in tIdx]) <= Si[int(i)]

def addFirstLastTaskConstraints(prob, numTasks, numDays, xihtVariables, xhitVariables):
    for t in range(numDays):
        xihtList = []
        xhitList = []
        for i in range(numTasks):
            xihtList.append(xihtVariables[i][t])
            xhitList.append(xhitVariables[i][t])
        prob += pulp.lpSum(xihtList) == pulp.lpSum(xhitList)
        prob += pulp.lpSum(xihtList) <= 1
        prob += pulp.lpSum(xhitList) <= 1

    def execute_algorithm(self, input_data_set):

        hosts = input_data_set.Hosts
        vms = input_data_set.VirtualMachines
        alert = input_data_set.Alert

        prob = pulp.LpProblem("test",pulp.LpMinimize)

        #if 1 host is enabled
        u = []

        x = []

        for host in hosts:

            ulp = pulp.LpVariable("used: %s" % host.Hostname,0,1,'Integer')

            u.append(ulp)

            innerX = []
            x.append(innerX)

            n = []
            m = []
            c = []

            for vm in vms:
                values = vm.getMetrics(host) #todo optimization

                lpVar = pulp.LpVariable("host %s contains %s" % (vm.InstanceName, host.Hostname),0,1,'Integer')
                innerX.append(lpVar)

                prob += ulp >= lpVar

                c.append([lpVar,values["C"]])
                n.append([lpVar,values["N"]])
                m.append([lpVar,values["M"]])

            prob += (pulp.LpAffineExpression(c) <= 1)
            prob += (pulp.LpAffineExpression(n) <= 1)
            prob += (pulp.LpAffineExpression(m) <= 1)


        #every vm on only one host
        for i in range(len(vms)):

            lps = [item[i] for item in x]
            prob += (pulp.lpSum(lps) == 1)


        prob += pulp.lpSum(u)
        GLPK().solve(prob)

        for v in prob.variables():
            print v.name, "=", v.varValue

def portfolio_lp(
    weights, latest_prices, min_allocation=0.01, total_portfolio_value=10000
):
    """
    For a long only portfolio, convert the continuous weights to a discrete allocation
    using Mixed Integer Linear Programming. This can be thought of as a clever way to round
    the continuous weights to an integer number of shares

    :param weights: continuous weights generated from the ``efficient_frontier`` module
    :type weights: dict
    :param latest_prices: the most recent price for each asset
    :type latest_prices: pd.Series or dict
    :param min_allocation: any weights less than this number are considered negligible,
                           defaults to 0.01
    :type min_allocation: float, optional
    :param total_portfolio_value: the desired total value of the portfolio, defaults to 10000
    :type total_portfolio_value: int/float, optional
    :raises TypeError: if ``weights`` is not a dict
    :raises TypeError: if ``latest_prices`` isn't a series
    :raises ValueError: if not ``0 < min_allocation < 0.3``
    :return: the number of shares of each ticker that should be purchased, along with the amount
             of funds leftover.
    :rtype: (dict, float)
    """
    import pulp

    if not isinstance(weights, dict):
        raise TypeError("weights should be a dictionary of {ticker: weight}")
    if not isinstance(latest_prices, (pd.Series, dict)):
        raise TypeError("latest_prices should be a pd.Series")
    if min_allocation > 0.3:
        raise ValueError("min_allocation should be a small float")
    if total_portfolio_value <= 0:
        raise ValueError("total_portfolio_value must be greater than zero")

    m = pulp.LpProblem("PfAlloc", pulp.LpMinimize)
    vals = {}
    realvals = {}
    etas = {}
    abss = {}
    remaining = pulp.LpVariable("remaining", 0)
    for k, w in weights.items():
        if w < min_allocation:
            continue
        vals[k] = pulp.LpVariable("x_%s" % k, 0, cat="Integer")
        realvals[k] = latest_prices[k] * vals[k]
        etas[k] = w * total_portfolio_value - realvals[k]
        abss[k] = pulp.LpVariable("u_%s" % k, 0)
        m += etas[k] <= abss[k]
        m += -etas[k] <= abss[k]
    m += remaining == total_portfolio_value - pulp.lpSum(realvals.values())
    m += pulp.lpSum(abss.values()) + remaining
    m.solve()
    results = {k: val.varValue for k, val in vals.items()}
    return results, remaining.varValue

    def add_stoichiometric_constraints(self, weights, S, compounds, reactions,
                                       net_reaction):
        """
            S is a NCxNR matrix where the rows represent compounds and the columns represent reactions.
            b is the linear constraint vector, and is NCx1
            
            compounds is the list of compounds from KEGG (NC long)
            reactions is a list of pairs of RID and direction (NR long)
        """
        self.weights = weights
        self.S = S
        self.compounds = compounds
        self.reactions = reactions
        
        assert S.shape[0] == len(self.compounds)
        assert S.shape[1] == len(self.reactions)

        self.cids = [c.cid for c in self.compounds]

        self.physiological_conc = np.matrix(np.ones((1, len(self.compounds)))) * default_c_mid
        if 1 in self.cids:
            self.physiological_conc[0, self.cids.index(1)] = 1 # [H2O] must be set to 1 in any case
        
        # reaction fluxes are the continuous variables
        self.flux_vars = pulp.LpVariable.dicts("Flux",
                                               [r.name for r in self.reactions],
                                               lowBound=0,
                                               upBound=self.flux_upper_bound,
                                               cat=pulp.LpContinuous)

        # add a linear constraint on the fluxes for each compound (mass balance)
        for c, compound in enumerate(self.compounds):
            reactions_with_c = list(np.nonzero(self.S[c, :])[1].flat)
            if len(reactions_with_c) == 0:
                continue
            
            # If the compound participates in the net reaction, force the mass
            # balance to be -1 times the desired net reaction coefficient
            if compound.cid in net_reaction.sparse:
                mass_balance = net_reaction.sparse[compound.cid]
            else:
                mass_balance = 0.0
                
            # Sum of all fluxes involving this compounds times the stoichiometry
            # of their reactions
            cid_sum = pulp.lpSum([self.S[c, r] * self.flux_vars[self.reactions[r].name]
                                  for r in reactions_with_c])
            
            self.prob.addConstraint(cid_sum == mass_balance,
                                    "C%05d_mass_balance" % compound.cid)
        
        obj = pulp.lpSum([self.flux_vars[self.reactions[r].name]*weight
                          for (r, weight) in self.weights])
        self.prob.setObjective(obj)

def opt_with_solver(node_ind, order_dict, travel_t, stay_t, pick_t, require_t, num,
                    mini_start, initial=None, load_check=True):
    # order_dict = {ord:(ori_id, dest_id)}, ord in order_set, ori_id, dest_id in node_ind
    # node_ind = range(n), travel_t = {(i,j): travel time between i and j}
    # initial = (given start id, given load)
    big_m = 10000
    eps = 1.0/10**7
    n = len(node_ind)
    inter_n = [(i, j) for i in node_ind for j in node_ind if j != i]
    off_time = lp.LpVariable('The route-off time')
    p = lp.LpVariable.dicts('Punish cost', node_ind, lowBound=0)
    x = lp.LpVariable.dicts('Route variable', inter_n, cat='Binary')
    o = lp.LpVariable.dicts('Start-point flag', node_ind, cat='Binary')
    d = lp.LpVariable.dicts('End-point flag', node_ind, cat='Binary')
    a = lp.LpVariable.dicts('Arrival time', node_ind)
    l = lp.LpVariable.dicts('Leave time', node_ind)
    t = lp.LpVariable.dicts('Order count', node_ind, 0, n-1+eps)
    if load_check:
        load = lp.LpVariable.dicts('Arrival load', node_ind, 0, MAX_LOADS+eps)
    else:
        load = lp.LpVariable.dicts('Arrival load', node_ind, 0)
    prob = lp.LpProblem('Optimize a route', lp.LpMinimize)
    # Objective
    prob += off_time + lp.lpSum(p[i] for i in node_ind)
    # Constraints
    for j in node_ind:
        prob += lp.lpSum(x[(i, j)] for i in node_ind if i != j) == 1 - o[j]
    for i in node_ind:
        prob += lp.lpSum(x[(i, j)] for j in node_ind if j != i) == 1 - d[i]
    prob += lp.lpSum(o[i] for i in node_ind) == 1
    prob += lp.lpSum(d[i] for i in node_ind) == 1
    if not (initial is None):
        prob += o[initial[0]] == 1
        prob += load[initial[0]] == initial[1]
    for order in order_dict:
        od = order_dict[order]
        prob += a[od[1]] >= l[od[0]]
    for i in node_ind:
        prob += off_time >= l[i]
        prob += l[i] >= a[i] + stay_t[i]
        prob += l[i] >= pick_t[i]
        prob += a[i] >= mini_start[i]
        prob += p[i] >= PUNISH_CO*(a[i] - require_t[i])
    for i, j in inter_n:
        prob += a[j] >= l[i] + travel_t[(i, j)] + big_m*(x[(i, j)] - 1)
        prob += t[j] >= t[i] + 1 + big_m*(x[(i, j)] - 1)
        prob += load[j] >= load[i] + num[i] + big_m*(x[(i, j)] - 1)
    prob.solve(lp.CPLEX(msg=1, timelimit=CPLEX_TIME_LIMIT, options=['set logfile cplex/cplex%d.log' % os.getpid()]))
    # set threads 100
    if lp.LpStatus[prob.status] != 'Infeasible':
        sol_list = [int(round(t[i].varValue)) for i in node_ind]
        return lp.value(prob.objective), sol_list, lp.LpStatus[prob.status]
    else:
        return None, None, lp.LpStatus[prob.status]