Python Model.getVars方法代码示例

# 需要导入模块: from gurobipy import Model [as 别名]
# 或者: from gurobipy.Model import getVars [as 别名]
def solve_lp_knapsack_gurobi(scores, costs, budget):
    from gurobipy import Model, LinExpr, GRB

    n = len(scores)

    # Create a new model.
    m = Model("lp_knapsack")

    # Create variables.
    for i in range(n):
        m.addVar(lb=0.0, ub=1.0)
    m.update()
    vars = m.getVars()

    # Set objective.
    obj = LinExpr()
    for i in range(n):
        obj += scores[i] * vars[i]
    m.setObjective(obj, GRB.MAXIMIZE)

    # Add constraint.
    expr = LinExpr()
    for i in range(n):
        expr += costs[i] * vars[i]
    m.addConstr(expr, GRB.LESS_EQUAL, budget)

    # Optimize.
    m.optimize()
    assert m.status == GRB.OPTIMAL
    x = np.zeros(n)
    for i in range(n):
        x[i] = vars[i].x

    return x

# 需要导入模块: from gurobipy import Model [as 别名]
# 或者: from gurobipy.Model import getVars [as 别名]
def tsp_gurobi(edges):
    """
    Modeled using GUROBI python example.
    """
    from gurobipy import Model, GRB, quicksum

    edges = populate_edge_weights(edges)
    incoming, outgoing, nodes = node_to_edge(edges)
    idx = dict((n, i) for i, n in enumerate(nodes))
    nedges = len(edges)
    n = len(nodes)

    m = Model()

    step = lambda x: "u_{0}".format(x)
    # Create variables
    vars = {}
    for i, (a, b, w) in enumerate(edges):
        vars[i] = m.addVar(obj=w, vtype=GRB.BINARY, name=str(i))
    for u in nodes[1:]:
        u = step(u)
        vars[u] = m.addVar(obj=0, vtype=GRB.INTEGER, name=u)
    m.update()

    # Bounds for step variables
    for u in nodes[1:]:
        u = step(u)
        vars[u].lb = 1
        vars[u].ub = n - 1

    # Add degree constraint
    for v in nodes:
        incoming_edges = incoming[v]
        outgoing_edges = outgoing[v]
        m.addConstr(quicksum(vars[x] for x in incoming_edges) == 1)
        m.addConstr(quicksum(vars[x] for x in outgoing_edges) == 1)

    # Subtour elimination
    edge_store = dict(((idx[a], idx[b]), i) for i, (a, b, w) in enumerate(edges))

    # Given a list of edges, finds the shortest subtour
    def subtour(s_edges):
        visited = [False] * n
        cycles = []
        lengths = []
        selected = [[] for i in range(n)]
        for x, y in s_edges:
            selected[x].append(y)
        while True:
            current = visited.index(False)
            thiscycle = [current]
            while True:
                visited[current] = True
                neighbors = [x for x in selected[current] if not visited[x]]
                if len(neighbors) == 0:
                    break
                current = neighbors[0]
                thiscycle.append(current)
            cycles.append(thiscycle)
            lengths.append(len(thiscycle))
            if sum(lengths) == n:
                break
        return cycles[lengths.index(min(lengths))]

    def subtourelim(model, where):
        if where != GRB.callback.MIPSOL:
            return
        selected = []
        # make a list of edges selected in the solution
        sol = model.cbGetSolution([model._vars[i] for i in range(nedges)])
        selected = [edges[i] for i, x in enumerate(sol) if x > .5]
        selected = [(idx[a], idx[b]) for a, b, w in selected]
        # find the shortest cycle in the selected edge list
        tour = subtour(selected)
        if len(tour) == n:
            return
        # add a subtour elimination constraint
        c = tour
        incident = [edge_store[a, b] for a, b in pairwise(c + [c[0]])]
        model.cbLazy(quicksum(model._vars[x] for x in incident) <= len(tour) - 1)

    m.update()

    m._vars = vars
    m.params.LazyConstraints = 1
    m.optimize(subtourelim)

    selected = [v.varName for v in m.getVars() if v.x > .5]
    selected = [int(x) for x in selected if x[:2] != "u_"]
    results = sorted(x for i, x in enumerate(edges) if i in selected) \
                    if selected else None
    return results

# 需要导入模块: from gurobipy import Model [as 别名]
# 或者: from gurobipy.Model import getVars [as 别名]

#.........这里部分代码省略.........
                        "Flow3[%s,%s,%s]" % (i, s, d),
                    )
        self.model.update()

    def Optimize(self, MipGap=None, TimeLimit=None, LogLevel=None):
        """ Optimize the defined  model.
    
        Parameters
        ----------
        MipGap : desired gap
        TimeLimit : time limit
        LogLevel: log level 1 for printing all optimal variables and None otherwise
        
        Returns
        -------
        BackupCapacity: The total capacity assigned per backup link
        BackupRoutes: The set of selected backup links 
           A tuple list with all paths for edge (s,d) that uses (i,j).
    
        """
        self.model.write("bpbackup.lp")

        if MipGap != None:
            self.model.params.MIPGap = MipGap
        if TimeLimit != None:
            self.model.params.timeLimit = TimeLimit
        # Compute optimal solution
        self.model.optimize()

        # Print solution
        if self.model.status == GRB.Status.OPTIMAL:

            if LogLevel == 1:
                for v in self.model.getVars():
                    print("%s %g" % (v.varName, v.x))

            self.BackupCapacitySolution = self.model.getAttr("x", self.BackupCapacity)
            self.BackupRoutesSolution = self.model.getAttr("x", self.bBackupLink)

            self.BackupLinksSolution = {}
            self.HatBackupCapacity = {}
            for link in self.BackupCapacitySolution:
                if self.BackupCapacitySolution[link] < 1 and self.BackupCapacitySolution[link] > 0.001:
                    self.HatBackupCapacity[link] = math.ceil(self.BackupCapacitySolution[link])
                else:
                    self.HatBackupCapacity[link] = math.floor(self.BackupCapacitySolution[link])
                if self.HatBackupCapacity[link] > 0:
                    if len(self.BackupLinksSolution) == 0:
                        self.BackupLinksSolution = [link]
                    else:
                        self.BackupLinksSolution = self.BackupLinksSolution + [link]
        else:
            print("Optimal value not found!\n")
            self.BackupCapacitySolution = []
            self.BackupRoutesSolution = {}
            self.BackupLinksSolution = {}

        return self.BackupCapacitySolution, self.BackupRoutesSolution, self.BackupLinksSolution, self.HatBackupCapacity

    def SaveBakupNetwork(self, file_name):
        """Save the optimal backup network to the file ``file_name``."""

        data = {
            "links": [i for i in self.BackupCapacitySolution],
            "capacities": [self.BackupCapacitySolution[i] for i in self.BackupCapacitySolution],
            "routes": [i for i in self.BackupRoutesSolution],

# 需要导入模块: from gurobipy import Model [as 别名]
# 或者: from gurobipy.Model import getVars [as 别名]

#.........这里部分代码省略.........
        if quadratic_component is not None:
            if not hasattr(quadratic_component, 'todok'):
                raise Exception('quadratic component must be a scipy.sparse type array')

            quadratic_objective = QuadExpr()
            for (index_0, index_1), the_value in quadratic_component.todok().items():
                quadratic_objective.addTerms(the_value,
                                       variable_list[index_0],
                                       variable_list[index_1])
            lp.setObjective(quadratic_objective, sense=objective_sense)
        #Constraints are based on mass balance
        #Construct the lin expression lists and then add
        #TODO: Speed this up as it takes about .18 seconds
        #HERE
        for the_metabolite in cobra_model.metabolites:
            constraint_coefficients = []
            constraint_variables = []
            for the_reaction in the_metabolite._reaction:
                constraint_coefficients.append(the_reaction._metabolites[the_metabolite])
                constraint_variables.append(reaction_to_variable[the_reaction])
            #Add the metabolite to the problem
            lp.addConstr(LinExpr(constraint_coefficients, constraint_variables),
                         sense_dict[the_metabolite._constraint_sense.upper()],
                         the_metabolite._bound,
                         the_metabolite.id)
    else:
        #When reusing the basis only assume that the objective coefficients or bounds can change
        if copy_problem:
            lp = the_problem.copy()
        else:
            lp = the_problem
        if not reuse_basis:
            lp.reset()
        for the_variable, the_reaction in zip(lp.getVars(),
                                              cobra_model.reactions):
            the_variable.lb = float(the_reaction.lower_bound)
            the_variable.ub = float(the_reaction.upper_bound)
            the_variable.obj = float(objective_sense*the_reaction.objective_coefficient)

    
    if the_problem == 'setup':
        return lp
    if print_solver_time:
        start_time = time()
    lp.update()
    lp.setParam("FeasibilityTol", tolerance_feasibility)
    lp.setParam("OptimalityTol", tolerance_optimality) 
    if tolerance_barrier:
        lp.setParam("BarConvTol", tolerance_barrier)

    if quad_precision:
            lp.setParam("Quad", 1)
    lp.setParam("Method", lp_method)

    #Different methods to try if lp_method fails
    the_methods = [0, 2, 1]
    if lp_method in the_methods:
        the_methods.remove(lp_method)
    if not isinstance(the_problem, Model):
        lp.optimize()
        if lp.status in status_dict:
            status = status_dict[lp.status]
        else:
            status = 'failed'
        if status != 'optimal':
            #Try to find a solution using a different method

# 需要导入模块: from gurobipy import Model [as 别名]
# 或者: from gurobipy.Model import getVars [as 别名]
    def generateInstance(self):
        # Check that variables have been intialized correctly
        if (not self.planes) or (self.n == 0) or (self.T == 0):
            logging.error('Store is not initialized correctly. Check for completeness of input-file.')
            return None

        model = Model("PlaneSolver")

        """
            r: Earliest landing possibility
            b: Best landing
            d: latest landing
            g: penalty for each time earlier from best
            h: penalty for each time later from best
            s_ij: if i lands before j: s needs to pass in time

            x_i: Landing time for i
            p_i: positive divertion from bi
            n_i: negative divertion from bi

        """

        bigM = max([max(x["s"]) for x in self.planes])
        logging.debug("Using bigM punishment: " + str(bigM))

        # generate Variables
        x = {}  # Landing time
        p = {}  # Positive divertion from optimal landing time
        n = {}  # Negative divertion from optimal landing time
        s = {}  # 1 iff plane i lands before j + buffer is invalid
        s1 = {}  # for these, see constraints
        s2 = {}
        h1 = {}
        h2 = {}
        for i in range(len(self.planes)):
            x[i] = model.addVar(name="x_%d" % (i+1), vtype="i", lb=self.planes[i]["r"], ub=self.planes[i]["d"])
            p[i] = model.addVar(name="p_%d" % (i+1), vtype="i", obj=self.planes[i]["h"], lb=0, ub=self.planes[i]["d"] - self.planes[i]["b"])
            n[i] = model.addVar(name="n_%d" % (i+1), vtype="i", obj=self.planes[i]["g"], lb=0, ub=self.planes[i]["b"] - self.planes[i]["r"])

            for j in range(len(self.planes)):
                s[i, j] = model.addVar(name="s_%d_%d" % (i+1, j+1), vtype="b", obj=bigM, lb=0, ub=1)
                s1[i, j] = model.addVar(name="s1_%d_%d" % (i+1, j+1), vtype="b", obj=0, lb=0, ub=1)
                s2[i, j] = model.addVar(name="s2_%d_%d" % (i+1, j+1), vtype="b", obj=0, lb=0, ub=1)
                h1[i, j] = model.addVar(name="h1_%d_%d" % (i+1, j+1), vtype="i", obj=0, lb=0)
                h2[i, j] = model.addVar(name="h2_%d_%d" % (i+1, j+1), vtype="i", obj=0, lb=0)

        model.modelSense = GRB.MINIMIZE
        model.update()

        for y in model.getVars():
            if y.ub < 10000:
                logging.debug("%10s\t[%5d:%5d]\tobj %4d", y.varName, y.lb, y.ub, y.obj)

        # Constraints
        for i in range(len(self.planes)):
            logging.debug("{} == {} + {} - {}".format(x[i].varName, self.planes[i]["b"], p[i].varName, n[i].varName))
            model.addConstr(x[i] == self.planes[i]["b"] + p[i] - n[i])

            for j in range(len(self.planes)):
                if j == i:
                    continue
                # model.addConstr(x[j] - x[i] )
                # Force s = s1 AND s2
                model.addConstr(s[i, j] <= s1[i, j])
                model.addConstr(s[i, j] <= s2[i, j])
                model.addConstr(s[i, j] >= s1[i, j] + s2[i, j] - 1)

                # s2 == xj - xi > 0
                logging.debug("{}\t <= {}\t - {}\t - {}\t".format(s2[i, j].varName, x[j].varName, x[i].varName, h1[i, j].varName))
                model.addConstr(s2[i, j] <= x[j] - x[i] + h1[i, j])
                logging.debug("{}\t >= 1\t".format(h1[i, j].varName))
                model.addConstr(h1[i, j] >= 1)
                logging.debug("{}\t - {}\t + {}\t <= {}\t*{}\t".format(x[j].varName, x[i].varName, h1[i, j].varName, s2[i, j].varName, bigM))
                model.addConstr(x[j] - x[i] + h1[i, j] <= s2[i, j]*bigM)

                # s1 == xi + sij - xj > 0
                logging.debug("{}\t + {}\t - {}\t >= {}\t".format(x[i].varName, self.planes[i]["s"][j], x[j].varName, s1[i, j].varName))
                model.addConstr(s1[i, j] <= x[i] + self.planes[i]["s"][j] - x[j] + h2[i, j])
                logging.debug("{}\t + {}\t - {}\t <= {}\t*{}\t".format(x[i].varName, self.planes[i]["s"][j], x[j].varName, s1[i, j].varName, bigM))
                model.addConstr(x[i] + self.planes[i]["s"][j] - x[j] <= s1[i, j]*bigM)

        # solve it
        model.optimize()

        # Debugging printouts
        dbgStrs = {}
        for y in model.getVars():
            if y.varName[0:2] not in dbgStrs:
                dbgStrs[y.varName[0:2]] = []
            dbgStrs[y.varName[0:2]].append("%8s = %4s [%4s]" % (y.varName, int(y.x), int(y.x*y.obj) if y.obj else ""))
        for x in dbgStrs:
            dbgStrs[x].sort()

        while dbgStrs:
            printed = False
            keys = [x for x in dbgStrs.keys()]
            keys.sort()
            logStr = ""
            for x in keys:
                if dbgStrs[x]:
#.........这里部分代码省略.........

# 需要导入模块: from gurobipy import Model [as 别名]
# 或者: from gurobipy.Model import getVars [as 别名]

#.........这里部分代码省略.........
        return val + y * width

    T = max(endTimes)

    x = {}
    wait = {}
    for i, j, t in itertools.product(range(n), range(width * height), range(T)):
        x[i, j, t] = model.addVar(name="x_%d_%d_%d" % (i+1, j+1, t+1), vtype="b", obj=1)
        wait[i, j, t] = model.addVar(name="wait_%d_%d_%d" % (i+1, j+1, t+1), vtype="i", obj=-1)

    model.modelSense = GRB.MINIMIZE
    model.update()

    # Force startpositions
    for i, (sN, sT) in enumerate(zip(startNodes, startTimes)):
        for t in range(sT):
            for j in range(width * height):
                if j == toGrid(sN - 1) and t == sT - 1:
                    logging.debug("start: %s == 1", x[i, j, t].varName)
                    model.addConstr(x[i, j, t] == 1)
                else:
                    logging.debug("start: %s == 0", x[i, j, t].varName)
                    model.addConstr(x[i, j, t] == 0)

    # Force endpositions
    for i, (eN, eT) in enumerate(zip(targetNodes, endTimes)):
        j = toGrid(eN - 1)
        logging.debug("end: %s == 1", x[i, j, eT - 1].varName)
        model.addConstr(x[i, j, eT - 1] == 1)
        # Container vanishes after endTime
        for t in range(eT, T):
            logging.debug("end: %s == 0", x[i, j, t].varName)
            model.addConstr(x[i, j, t] == 0)

    # single container per node
    for j, t in itertools.product(range(width * height), range(T)):
        logging.debug("%s <= 1", [x[i, j, t].varName for i in range(n)])
        model.addConstr(quicksum(x[i, j, t] for i in range(n)) <= 1)

    # Force valid container movement
    for w, h in itertools.product(range(width), range(height)):
        vals = [toGrid(w, h)]
        if h >= 1:
            vals += [toGrid(w, h - 1)]
        if w >= 1:
            vals += [toGrid(w - 1, h)]
        if h+1 < height:
            vals += [toGrid(w, h + 1)]
        if w+1 < width:
            vals += [toGrid(w + 1, h)]

        for i, t in itertools.product(range(n), range(1, T)):
            if endTimes[i] > t and startTimes[i] <= t:
                logging.debug("sum(%s) >= %s", [x[i, j, t].varName for j in vals], x[i, toGrid(w, h), t - 1].varName)
                model.addConstr(quicksum(x[i, j, t] for j in vals) >= x[i, toGrid(w, h), t - 1])
            else:
                logging.debug("skipped(%s) >= %s", [x[i, j, t].varName for j in vals], x[i, toGrid(w, h), t - 1].varName)

        for i, t in itertools.product(range(n), range(1, T)):
            logging.debug("sum(%s) <= 1", [x[i, j, t].varName for j in vals])
            model.addConstr(quicksum(x[i, j, t] for j in vals) <= 1)

    # Force continous line through grid
    for i, t in itertools.product(range(n), range(0, T)):
        if endTimes[i] > t and startTimes[i] <= t + 1:
            logging.debug("sum(%s) == 1", [x[i, j, t].varName for j in range(width * height)])
            model.addConstr(quicksum(x[i, j, t] for j in range(width * height)) == 1)
        else:
            logging.debug("sum(%s) == 0", [x[i, j, t].varName for j in range(width * height)])
            model.addConstr(quicksum(x[i, j, t] for j in range(width * height)) == 0)

    # Prevent ships from passing over same link
    for t in range(1, T):
        for w, h in itertools.product(range(width - 1), range(height)):
            for i in range(n):
                for k in range(i+1, n):
                    model.addConstr(x[i, toGrid(w, h), t - 1] + x[i, toGrid(w + 1, h), t] + x[k, toGrid(w, h), t] + x[k, toGrid(w + 1, h), t - 1] <= 3)
                    model.addConstr(x[i, toGrid(w, h), t] + x[i, toGrid(w + 1, h), t - 1] + x[k, toGrid(w, h), t - 1] + x[k, toGrid(w + 1, h), t] <= 3)
        for w, h in itertools.product(range(width), range(height - 1)):
            for i in range(n):
                for k in range(i+1, n):
                    model.addConstr(x[i, toGrid(w, h), t - 1] + x[i, toGrid(w, h + 1), t] + x[k, toGrid(w, h), t] + x[k, toGrid(w, h + 1), t - 1] <= 3)
                    model.addConstr(x[i, toGrid(w, h), t] + x[i, toGrid(w, h + 1), t - 1] + x[k, toGrid(w, h), t - 1] + x[k, toGrid(w, h + 1), t] <= 3)

    # Allow free waiting
    for i, j, t in itertools.product(range(n), range(width * height), range(0, T)):
        if t < startTimes[i]:
            model.addConstr(x[i, j, t] == wait[i, j, t])
        else:
            model.addConstr(x[i, j, t - 1] + x[i, j, t] - 1 <= wait[i, j, t])
            model.addConstr(x[i, j, t - 1] >= wait[i, j, t])
            model.addConstr(x[i, j, t] >= wait[i, j, t])

    model.optimize()

    for y in model.getVars():
        if y.x:
            logging.warning("%s = %d", y.varName, y.x)

    return model

# 需要导入模块: from gurobipy import Model [as 别名]
# 或者: from gurobipy.Model import getVars [as 别名]

#.........这里部分代码省略.........
        
    def __init__(self,imp_samp,nodes,links,capacity,epsilon,N,backup_link,link_capacity):
        '''
        Constructor
        '''
        self.__links = links
        self.__nodes = nodes
        self.__capacity = capacity
        self.__epsilon = epsilon
        self.__N = N
        self.__loadModel(imp_samp,backup_link,link_capacity)
                                 
    def __loadModel(self,imp_samp, backup_link,link_capacity):
                 
        # Create optimization model
        self.__model = Model('Backup')
     
        for i,j in self.__links:
            for k in range(self.__N):
                self.__z[k,i,j] = self.__model.addVar(lb=0,name='z[%s][%s][%s]' % (k,i,j))
        self.__model.update()
        
        for i,j in self.__links: 
            self.__z0[i,j] = self.__model.addVar(lb=-GRB.INFINITY,name='z0[%s][%s]' %(i,j))
        self.__model.update()
         
        for i,j in self.__links: 
            self.__q[i,j] = self.__model.addVar(lb=-GRB.INFINITY,name='q[%s][%s]' %(i,j))
        self.__model.update()
        
        self.__model.modelSense = GRB.MINIMIZE
        
        self.__model.setObjective(quicksum(self.__q[i,j] for i,j in self.__links))
        
        self.__model.update()
         
            
        #------------------------------------------------------------------------#
        #                    Constraints definition                              #
        #                                                                        #
        #                                                                        #
        #------------------------------------------------------------------------#
          
        # Buffer probability I
        for i,j in self.__links:
            self.__model.addConstr(self.__z0[i,j] + 1/(self.__N*self.__epsilon)*quicksum(self.__z[k,i,j]*imp_samp[k] for (k) in range(self.__N)) <= self.__q[i,j],'[CONST]Buffer_Prob_I[%s][%s]'%(i,j))
        self.__model.update()
         
        # Link capacity constraints
        for i,j in self.__links:
            for k in range(self.__N):
                self.__model.addConstr((quicksum(backup_link[i,j,s,d]*self.__capacity[k,s,d] for s,d in self.__links) - link_capacity[i,j] - self.__z0[i,j]) <= self.__z[k,i,j],'[CONST]Buffer_Prob_II[%s][%s][%s]' % (k,i,j))
        self.__model.update()
        
        # Link capacity constraints
        for i,j in self.__links:
            for k in range(self.__N):
                self.__model.addConstr(self.__z[k,i,j] >= 0,'[CONST]Buffer_Prob_III[%s][%s][%s]' % (k,i,j))
        self.__model.update()
         
    def optimize(self,MipGap, TimeLimit, LogLevel = None):
         
        self.__model.write('quantile.lp')
         
        if MipGap != None:
            self.__model.params.MIPGap = MipGap
        if TimeLimit != None:
            self.__model.params.timeLimit = TimeLimit
        # Compute optimal solution
        self.__model.optimize()
         
        # Print solution
        if self.__model.status == GRB.Status.OPTIMAL:
            #SuperQuantileSolution = self.__model.getAttr('x', self.__z0)
            SuperQuantileSolution = {}
            OptimalZnot = {}
            for i,j in self.__links:
                name='q[%s][%s]'%(i,j)
                v = self.__model.getVarByName(name)
                SuperQuantileSolution[i,j]=v.x
                name='z0[%s][%s]'%(i,j)
                v = self.__model.getVarByName(name)
                OptimalZnot[i,j]=v.x
            
            if LogLevel == 1:
                for v in self.__model.getVars():
                    print('%s %g' % (v.varName, v.x))
                                                
        else:
            print('Optimal value not found!\n')
            SuperQuantileSolution = {}
            OptimalZnot={}
                         
        return SuperQuantileSolution, OptimalZnot    
    
    def reset(self):
        '''
        Reset model solution
        '''
        self.__model.reset()

# 需要导入模块: from gurobipy import Model [as 别名]
# 或者: from gurobipy.Model import getVars [as 别名]
    def createModel(self):
        w = {}
        self.y = {}
        
        m = Model("Optimization Model")
        
        #part of 6a -- create w variables
        for n in range(1,self.nScenario+1):
            for k in range(1,self.numberOfFinancialAsstValues+1): #I'm not sure how the "paramDF" file will be structured--this is temporary
                for j in range(1,self.nOwners+1):
                    w[j, k, n] = m.addVar(vtype=GRB.CONTINUOUS, name="w_"+str(j)+"_"+str(k)+"_"+str(n))   

        #Constraint 6f
        for k in range(1,self.numberOfFinancialAsstValues+1):
            for j in range(1,self.nOwners+1):
                self.y[j, k] = m.addVar(vtype=GRB.BINARY, name="y_"+str(j)+"_"+str(k))

        m.update()
        
        #6a continued
        #for k in range(1,self.numberOfFinancialAsstValues+1):
        #6a updated
        for n in range(1,self.nScenario+1):
            m.addConstr(quicksum(w[1,k,n] for k in range(1,self.numberOfFinancialAsstValues+1)) == self.SecondStgValues[n-1]*quicksum(self.DecisionProb[n,1,k]*self.y[1, k] for k in range(1,self.numberOfFinancialAsstValues+1)), name = "6a_1_"+str(n))
        #        w[1,k,n] = self.SecondStgValues[n-1]
        
        #6b updated
        for r in range(2,self.nOwners+1):
            for n in range(1,self.nScenario+1):
                #if self.ProbDict["("+str(n)+", 1, 0, "+str(k)+")"] > 0:
        #        for self.ProbDict["("+str(n)+", 1, 0, "+str(k)+")"] > 0:
                    m.addConstr(quicksum(w[r-1,k,n] for k in range(1,self.numberOfFinancialAsstValues+1)) == quicksum(w[r,k,n]*(1/self.DecisionProb[n,r,k]) for k in range(1,self.numberOfFinancialAsstValues+1)), name = "6b_"+str(r)+"_"+str(n))
                #if self.ProbDict["("+str(n)+", 1, 1, "+str(k)+")"] > 0:
                #    m.addConstr(quicksum(w[r-1,k,n] for k in range(1,self.numberOfFinancialAsstValues+1)) == quicksum(w[r,k,n]*(1/self.ProbDict["("+str(n)+", 1, 1, "+str(k)+")"]) for k in range(1,self.numberOfFinancialAsstValues+1)), name = "6b_"+str(r)+"_"+str(n))

        #6b
        #for r in range(2, self.nOwners+1):
        #    for n in range(1,self.nScenario+1):
        #        m.addConstr(quicksum(quicksum(self.ProbDict["("+str(n)+", "+str(r-1)+", "+str(l)+", "+str(k)+")"]*w[r-1,k,n] for l in (0,1))
        #                             for k in range(1,self.numberOfFinancialAsstValues+1)) == quicksum(w[r, k, n] for k in range(1,self.numberOfFinancialAsstValues+1)), name = "6b_"+str(r)+"_"+str(n))
        
        #for n in range(1,self.nScenario+1):
        #    m.addConstr(quicksum(w[r, k, n] for k in range(1,self.numberOfFinancialAsstValues+1)), name = "6b_"+str(r)+"_"+str(n))
        
        
        #for k in range(1,self.numberOfFinancialAsstValues+1):
        #    for r in range(2, self.ownerNums+2):
        #        for n in range(1,self.nScenario+1):
        #            for l in (0,1):
        #                m.addConstr(quicksum(quicksum(self.ProbDict["("+str(n)+", "+str(r-1)+", "+str(l)+", "+str(k)+")"]*w[r-1,k,n]) == quicksum(w[r, k, n]),
        #                                     name = "6b_"+str(r)+"_"+str(k)+"_"+str(n)))
                
                
        #6c
        #Not sure if this is the proper formatting for this constraint
        #Is this the proper use of self.SecondStgValues?
        #for n in range(1,self.nScenario +1):
        #    for k in range(1,self.numberOfFinancialAsstValues+1): #I'm not sure how the "paramDF" file will be structured--this is temporary
        #        for r in range(1,self.ownerNums+2):
        #6c updated            
        for k in range(1,self.numberOfFinancialAsstValues+1):
            for r in range(1, self.nOwners+1):
                for n in range(1,self.nScenario+1):
                    m.addConstr(w[r, k, n] <= self.y[r, k]*self.SecondStgValues[n-1], name = "6c_"+str(r)+"_"+str(k)+"_"+str(n))
                    #print str(r)+"_"+str(k)+"_"+str(n)
                                
                                
                                #for r in range(1, self.ownerNums+2) for k in range(1,self.numberOfFinancialAsstValues+1) for n in range(1,self.nScenario+1))
        
        #Constraint 6d
        #the sum of the financial assistance offered to all landowners is less than or equal to the agency's budget
        #Where does C come from?
        m.addConstr(quicksum(quicksum(self.C_k[k-1]*self.y[j, k] for k in range(1,self.numberOfFinancialAsstValues+1))for j in range(1,self.nOwners+1)) <= self.Budget_param, name = "6d")
        
        #Constraint 6e
        for j in range(1,self.nOwners+1):
            m.addConstr(quicksum(self.y[j, k] for k in range(1,self.numberOfFinancialAsstValues+1)) == 1, name = "6e_"+str(j))
                        
        m.update()
                                              
        #set objective
        lastLandownerIndex = self.nOwners
        
        m.setObjective(quicksum(quicksum(w[lastLandownerIndex, k, n] for k in range(1,self.numberOfFinancialAsstValues+1)) for n in range(1,self.nScenario+1)), GRB.MINIMIZE)        

        m.update()
        
        m.optimize()
        
        if m.status == GRB.Status.OPTIMAL:
            print ('\nOBJECTIVE VALUE: %g' % m.objVal)
            
        for v in m.getVars():
            print('%s %g' % (v.varName, v.x))
       
        m.write('toy results.lp')
        
        return m

# 需要导入模块: from gurobipy import Model [as 别名]
# 或者: from gurobipy.Model import getVars [as 别名]
class HarmonyModel:

    def __init__(self, length):

        """Initialize model, note and chord variables, and constraints for notes
        and chord constraints."""

        self.model = Model("harmony")
        self.length = length

        config = load_source("config",
                join(abspath(dirname(__file__)), "config.py"))

        config.length = self.length

        config.model = self.model

        self.notes = Notes(config)
        config.notes = self.notes

        self.chords = Chords(config)


    def ensure_melody(self, melody):

        """Ensure that the four-part harmony has a complete or partial melody as
        its soprano line."""

        for t, note in enumerate(melody):

            if not note:
                continue

            constraint = self.notes[denote(note), 3, t] == 1

            self.model.addConstr(constraint,
                "melody_t" + str(t))

            self.model.update()


    def ensure_harmony(self, harmony):

        """Ensure that the four-part harmony uses certain chords (and the
        correct inversion) at each time step."""

        for t, chord in enumerate(harmony):

            if not chord:
                continue

            for c in self.chords.chord_notes.keys():

                if c == chord:
                    continue

                for d in self.chords.chord_doublings[c]:

                    constraint = self.chords[c, d, t] == 0

                    self.model.addConstr(constraint,
                        "harmony_c" + c + "_d" + str(d) + "_t" + str(t))

            self.model.update()


    def write(self):

        """Write out LP file of generated model (for debugging)."""

        self.model.write("model.lp")


    def solve(self):

        """Solve the model and provide the solution's notes (x variables) and
        chords (y variables)."""

        self.model.optimize()

        self.solution = []

        try:

            for var in self.model.getVars():
                if var.x == 1 and var.varName[0] in {"x", "y"}:
                    self.solution.append(var.varName)

            return self.solution

        except:

            return []
        

    def lilypond_export(self):

        """Return a Lilypond file to generate a human-readable score based on
        the IP solution."""

#.........这里部分代码省略.........

# 需要导入模块: from gurobipy import Model [as 别名]
# 或者: from gurobipy.Model import getVars [as 别名]
    def generateInstance(self):
        # Check that store has been intialized correctly
        if not all([x in self.store for x in ["Timehorizon", "h", "K", "a", "d", "s", "st"]]):
            logging.error('Store is not initialized correctly. Check for completeness of input-file.')
            return None

        model = Model("LotSolver")

        tHor = self.store["Timehorizon"]
        nPr = self.store["nProducts"]

        timeRange = range(1, tHor + 1)
        prodRange = range(1, nPr + 1)

        # Assume that production of products costs at least 1h, so a max of K products can be produced per period
        bigM = max(self.store["K"])

        # generate Variables
        # boolean production
        bx = {}
        # lager
        l = {}
        # production
        x = {}
        for t in timeRange:
            for p in prodRange:
                x[p, t] = model.addVar(name="x_%d_%d" % (p, t), vtype="i")
                bx[p, t] = model.addVar(name="bx_%d_%d" % (p, t), vtype="b")
        for t in range(0, tHor + 1):
            for p in prodRange:
                l[p, t] = model.addVar(name="l_%d_%d" % (p, t), vtype="i", obj=float(self.store["h"][p-1]))

        # switch costs
        s = {}
        for t in range(0, tHor+1):
            for p1 in prodRange:
                for p2 in prodRange:
                    # Switching to the same product does not cost anything - even if the file may say otherwise
                    objective = float(self.store["s"][p1-1][p2-1]) if p1 != p2 else 0
                    s[p1, p2, t] = model.addVar(name="s_%d_%d_%d" % (p1, p2, t), vtype="b", obj=objective)

        model.modelSense = GRB.MINIMIZE
        model.update()

        for y in model.getVars():
            logging.debug("%s obj %s", y.varName, y.obj)

        # Constraints

        # Initially only allow a single product
        logging.debug("%s == 1", [s[key].varName for key in s if key[2] == 0])
        model.addConstr(quicksum(s[key] for key in s if key[2] == 0) == 1, name="single_switch_" + str(t))

        # Only allow products in each period that actually has been switched to
        for t in timeRange:
            for p in prodRange:
                logging.debug("(%s) == %s", [s[key].varName for key in s if key[2] == t and (key[1] == p or key[0] == p)], bx[p, t].varName)
                model.addConstr(quicksum(s[key] for key in s if key[2] == t and (key[1] == p or key[0] == p)) >= bx[p,t], name="single_switch_" + str(t))

                # Force bx = 1 iff x != 0
                model.addConstr(x[p,t] >= bx[p,t])
                model.addConstr(x[p,t] <= bigM * bx[p,t])

        for t in timeRange:
            # Allow only a single switch each period
            logging.debug('Single switch constraint for ' + str([s[key].varName for key in s if key[2] == t]))
            model.addConstr(quicksum(s[key] for key in s if key[2] == t) == 1, name="single_switch_" + str(t))

            # Only allow connected switches between t-1 and t
            for p in prodRange:
                logging.debug('valid_switch for ' + str([s[key].varName for key in s if key[2] == (t-1) and key[1] == p]) + " and " + str([s[key].varName for key in s if key[2] == t and key[0] == p]))
                model.addConstr(quicksum(s[key] for key in s if key[2] == (t-1) and key[1] == p) == quicksum(s[key] for key in s if key[2] == t and key[0] == p), 
                    name="valid_switch_%d_%d" % (p, t))

        # Machine can't be occupied for more then K hours / period
        for t in timeRange:
            logging.debug("sum {} + sum {} <= {}".format([x[key].varName + "*" + str(self.store["a"][key[0]-1]) for key in x if key[1] == t],
                [s[key].varName + "*" + str(self.store["st"][key[0]-1][key[1]-1]) for key in s if key[2] == t], self.store["K"][t-1]))
            model.addConstr(quicksum(x[key]*self.store["a"][key[0]-1] for key in x if key[1] == t) + quicksum(s[key]*self.store["st"][key[0]-1][key[1]-1] for key in s if key[2] == t) <= self.store["K"][t-1])

        # Initial warehouse stock
        for p in prodRange:
            logging.debug("%s == %s", l[p, 0].varName, self.store["l"][p-1])
            model.addConstr(l[p, 0] == self.store["l"][p-1])
        # Update warehouse stock inbetween periods + enforce demand to be met
        for t in range(1, tHor+1):
            for p in prodRange:
                logging.debug("{} = {} + {} - {}".format(l[p, t].varName, l[p, t-1].varName, x[p, t].varName, self.store["d"][p-1][t-1]))
                model.addConstr(l[p, t] == l[p, t-1] + x[p, t] - self.store["d"][p-1][t-1])

        # solve it
        model.optimize()

        # Debugging printouts
        for y in model.getVars():
            if y.x >= 0.001:
                logging.debug("%s = %s cost %d", y.varName, y.x, y.x*y.obj)

        for t in timeRange:
            logging.debug("%s + %s", (["{}[{}] * {}".format(x[key].x, x[key].varName, self.store["a"][key[0]-1]) for key in x if key[1] == t]), ([str(s[key].varName) + "*" + str(self.store["st"][key[0]-1][key[1]-1]) for key in s if key[2] == t and s[key].x >= 0.001]))           
#.........这里部分代码省略.........