Python gurobipy.quicksum() Examples

def fit(array, convex=1):
    """Fit a smooth line to the given time-series data"""
    N = len(array)
    m = gurobipy.Model()
    fv = m.addVars(N)
    if convex == 1:
        m.addConstrs(fv[i] <= fv[i-1] for i in range(1,N))
        m.addConstrs(fv[i] + fv[i-2] >= 2*fv[i-1] for i in range(2,N))
    else:
        m.addConstrs(fv[i] >= fv[i-1] for i in range(1,N))
        m.addConstrs(fv[i] + fv[i-2] <= 2*fv[i-1] for i in range(2,N))
    m.setObjective(
        gurobipy.quicksum([fv[i] * fv[i] for i in range(N)])
        - 2 * gurobipy.LinExpr(array,fv.values())
    )
    m.Params.outputFlag = 0
    m.optimize()
    return [fv[i].X for i in range(N)] 

def _subtourelim(model, where):
  """  Callback - use lazy constraints to eliminate sub-tours """
  if where == GRB.callback.MIPSOL:
    # make a list of edges selected in the solution
    X = model.cbGetSolution(model._vars)
    n = int(sqrt(len(X)))
    selected = [(i,j) for i in xrange(n) for j in xrange(n) if X[(i,j)]>0.5]

    # find the shortest cycle in the selected edge list
    tour = _subtour(selected,n)
    if len(tour) < n:
      # add a subtour elimination constraint
      expr = quicksum(model._vars[tour[i], tour[j]]
          for i in xrange(len(tour))
              for j in xrange(i+1, len(tour)))
      model.cbLazy(expr <= len(tour)-1)

# closures that include n 

def _add_set_covering_constraints(m, N, K, X_j, X_j_keys, X_j_node_sets, X_j_forbidden_combos):
    # c1, the convering constraint, each node mus be served exactly by 
    #  one route
    c1_constrs = []
    for i in xrange(1,N):
        active_vars_sum = quicksum([X_j[j]
                                   for j, ns in enumerate(X_j_node_sets)
                                   if i in ns])
        c1c = m.addConstr( active_vars_sum == 1, "c1_n%d"%i )
        c1_constrs.append(c1c);
        
    # c2, the forbid constraints that do not allow some petal combinations
    c2_constrs = []
    for i, fc in enumerate(X_j_forbidden_combos):
        # sum_{i \in S } x_i - sum_{j \in \not{S} } x_j < |S| 
        combo_sum = quicksum([X_j[j] if (j in fc) else -X_j[j]
                                    for j in X_j_keys])
        c2c = m.addConstr( combo_sum  <= len(fc)-1, "c2_n%d"%i )
        c2_constrs.append(c2c);
    
    # c3, number of vehicles constraint
    c3_constrs = []
    if K:
        active_vars_sum = quicksum(X_j)
        c3c = m.addConstr( active_vars_sum <= K, "c3" )
        c3_constrs.append(c3c)
        
    return c1_constrs, c2_constrs, c3_constrs 

def _set_objective_function(self):
        # Similar to constraints, saving the costs expressions as attributes
        # can give you the chance to retrieve their values at the end of the optimization
        self.total_holding_cost = self.input_params['holding_cost'] * grb.quicksum(self.inventory_variables.values())

        self.total_production_cost = grb.quicksum(row['production_cost'] * self.production_variables[index]
                                                  for index, row in self.input_data.iterrows())

        objective = self.total_holding_cost + self.total_production_cost
        self.model.setObjective(objective, grb.GRB.MINIMIZE)

    # ================== Optimization ================== 

def regularize(self, center, norm, a, b, i):
        """Regularize a stochastic model.

        Parameters
        ----------

        center: array-like
            The regularization center with length n_states.

        norm: 'L1'/'L2'
            The norm to use for regularization.

        a,b,i: float,float,integer (a>0, 0<b<1, i>0)
            The coefficient of the regularization term is a*b^{i}, where i is
            the index of iteration.
        """
        self.rgl = self._model.addVar(
            lb=0,
            obj=self.modelsense*b**i,
            name='rgl'
        )
        if norm == 'L1':
            self.rgl_constr = self._model.addConstrs(
                (self.rgl >= a*(self.states[i] - center[i])
                for i in range(self.n_states)),
                name = 'rgl'
            ).values()
        elif norm == 'L2':
            self.rgl_constr = [self._model.addQConstr(
                self.rgl -
                a*gurobipy.QuadExpr(
                    gurobipy.quicksum([
                        self.states[i] * self.states[i]
                        - self.states[i] * 2 * center[i]
                        + center[i] * center[i]
                        for i in range(self.n_states)
                    ])
                )
                >=0,
                name = 'rgl'
            )]
        self._model.update() 

def _binarize(self, precision, n_binaries, transition=0):
        # Binarize StochasticModel. StochasticModel at transition stage keeps
        # states in original space while binarzing local_copies
        self.n_states_original_space = self.n_states
        self.local_copies_original_space = self.local_copies
        self.states_original_space = self.states
        if transition == 0:
            self.states = []
            self.n_states = 0
        self.local_copies = []
        for i, (x, y) in enumerate(
            zip(self.states_original_space, self.local_copies_original_space)
        ):
            if transition == 0:
                states = self.addVars(
                    n_binaries[i], vtype=gurobipy.GRB.BINARY, name=x.varName
                ).values()
            local_copies = self.addVars(
                n_binaries[i], vtype=gurobipy.GRB.BINARY, name=y.varName
            ).values()
            self.update()
            if transition == 0:
                temp1 = gurobipy.quicksum(
                    pow(2, i) * states[i] for i in range(n_binaries[i])
                )
            temp2 = gurobipy.quicksum(
                pow(2, i) * local_copies[i] for i in range(n_binaries[i])
            )
            # Assume bounds are the same over time!
            if x.vtype not in ["I","B"]:
                if transition == 0:
                    self.addConstr(
                        temp1 == precision * (x - x.lb),
                        name="binarize_states_{}".format(i),
                    )
                self.addConstr(
                    temp2 == precision * (y - y.lb),
                    name="binarize_local_copies_{}".format(i),
                )
            else:
                x.lb = math.ceil(x.lb)
                y.lb = math.ceil(y.lb)
                if transition == 0:
                    self.addConstr(
                        temp1 == x - x.lb,
                        name="binarize_states_{}".format(i),
                    )
                self.addConstr(
                    temp2 == y - y.lb,
                    name="binarize_local_copies_{}".format(i),
                )
            if transition == 0:
                self.states += states
                self.n_states += n_binaries[i]
            self.local_copies += local_copies 

def gurobi_tsp(distance_matrix):
    """
    Solves tsp problem.
    :param distance_matrix: symmetric matrix of distances, where the i,j element is the distance between object i and j
    :return: matrix containing {0, 1}, 1 for each transition that is included in the tsp solution
    """
    n = len(distance_matrix)
    m = Model()
    m.setParam("OutputFlag", False)
    m.setParam("Threads", 1)

    # Create variables
    vars = {}
    for i in range(n):
        for j in range(i + 1):
            vars[i, j] = m.addVar(
                obj=0.0 if i == j else distance_matrix[i][j], vtype=GRB.BINARY, name="e" + str(i) + "_" + str(j)
            )
            vars[j, i] = vars[i, j]
        m.update()

    # Add degree-2 constraint, and forbid loops
    for i in range(n):
        m.addConstr(quicksum(vars[i, j] for j in range(n)) == 2)
        vars[i, i].ub = 0
    m.update()

    # Optimize model
    m._vars = vars
    m.params.LazyConstraints = 1

    def subtour_fn(model, where):
        return subtourelim(n, model, where)

    m.optimize(subtour_fn)
    solution = m.getAttr("x", vars)
    selected = [(i, j) for i in range(n) for j in range(n) if solution[i, j] > 0.5]
    result = np.zeros_like(distance_matrix)
    for (i, j) in selected:
        result[i][j] = 1

    return result 

def create_model(instance):
  """
  Creates a simple MIP model for a set cover instance.

  Creates one binary decision variable s_i for each set. The objective
  function is simply the sum over the c_i * s_i. The constraints
  are that each item is captured by at least one set that is taken.

  Args:
    instance: The set cover instance as created by read().

  Returns:
    A pair of the Gurobi MIP model and the mapping from the sets
    in the instance to the corresponding Gurobi variables.
  """
  name, nitems, sets = instance
  model = grb.Model(name)

  # One variable for each set. Also remember which sets cover each item.
  covered_by = [[] for i in range(nitems)]
  vars = []
  for i, set in enumerate(sets):
    cost, covers = set
    vars.append(model.addVar(obj=cost, vtype=grb.GRB.BINARY, name="s_{0}".format(i)))

    for item in covers:
      covered_by[item].append(vars[i])
  model.update()

  # Constraint: Each item covered at least once.
  for item in range(nitems):
    model.addConstr(grb.quicksum(covered_by[item]) >= 1)

  # We want to minimize. Objective coefficients already fixed during variable creation.
  model.setAttr("ModelSense", grb.GRB.MINIMIZE)

  # Tuning parameters derived from sc_330_0
  model.read("mip.prm")

  model.setParam("Threads", 3)
  model.setParam("MIPGap", 0.001)  # 0.1% usually suffices

  return model, vars 

def create_model(instance):
  """
  Creates a simple MIP model for a set cover instance.

  Creates one binary decision variable s_i for each set. The objective
  function is simply the sum over the c_i * s_i. The constraints
  are that each item is captured by at least one set that is taken.

  Args:
    instance: The set cover instance as created by read().

  Returns:
    A pair of the Gurobi MIP model and the mapping from the sets
    in the instance to the corresponding Gurobi variables.
  """
  name, nitems, sets = instance
  model = grb.Model(name)

  # One variable for each set. Also remember which sets cover each item.
  covered_by = [[] for i in range(nitems)]
  vars = []
  for i, set in enumerate(sets):
    cost, covers = set
    vars.append(model.addVar(obj=cost, vtype=grb.GRB.BINARY, name="s_{0}".format(i)))

    for item in covers:
      covered_by[item].append(vars[i])
  model.update()

  # Constraint: Each item covered at least once.
  for item in range(nitems):
    model.addConstr(grb.quicksum(covered_by[item]) >= 1)

  # We want to minimize. Objective coefficients already fixed during variable creation.
  model.setAttr("ModelSense", grb.GRB.MINIMIZE)

  # Tuning parameters derived from sc_330_0
  model.read("mip.prm")

  model.setParam("Threads", 3)
  model.setParam("MIPGap", 0.001)  # 0.1% usually suffices

  return model, vars