Python Model.getConstrs方法代码示例

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

#.........这里部分代码省略.........
                         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
            lp.setParam("MarkowitzTol", 1e-2)
            for lp_method in the_methods:
                lp.setParam("Method", lp_method)
                lp.optimize()
                if status_dict[lp.status] == 'optimal':
                    break
    else:
        lp.setParam("TimeLimit", 0.6)
        lp.optimize()
        lp.setParam("TimeLimit", "default")
        if lp.status in status_dict:
            status = status_dict[lp.status]
        else:
            status = 'failed'
        if status != 'optimal':
            lp.setParam("MarkowitzTol", 1e-2)
            #Try to find a solution using a different method
            for lp_method in the_methods:
                lp.setParam("Method", lp_method)
                lp.optimize()
                if status_dict[lp.status] == 'optimal':
                    break
                
            if status_dict[lp.status] != 'optimal':
                lp = optimize_gurobi(cobra_model, new_objective=new_objective, objective_sense=objective_sense,
                                     min_norm=min_norm, the_problem=None, 
                                     print_solver_time=print_solver_time)['the_problem']


    if print_solver_time:
        print 'optimize time: %f'%(time() - start_time)
    x_dict = {}
    y_dict = {}
    y = None
    if lp.status in status_dict:
        status = status_dict[lp.status]
    else:
        status = 'failed'
    if status == 'optimal':
        objective_value = objective_sense*lp.ObjVal
        [x_dict.update({v.VarName: v.X}) for v in lp.getVars()]
        x = array([x_dict[v.id] for v in cobra_model.reactions])
        if lp.isMIP:
            y = y_dict = None #MIP's don't have duals
        else:
            [y_dict.update({c.ConstrName: c.Pi})
             for c in lp.getConstrs()]
            y = array([y_dict[v.id] for v in cobra_model.metabolites])
    else:
        y = y_dict = x = x_dict = None
        objective_value = None
        if error_reporting:
            print 'gurobi failed: %s'%lp.status  
    the_solution = Solution(objective_value, x=x, x_dict=x_dict,
                            y=y, y_dict=y_dict,
                            status=status)
    solution = {'the_problem': lp, 'the_solution': the_solution}
    return solution

# 需要导入模块: from gurobipy import Model [as 别名]
# 或者: from gurobipy.Model import getConstrs [as 别名]
def _sensitivity_for_constraints(AT, j, project, y_, project_activity, M):
    global _round, _pa_dataset
    m = Model("SingleProject_%d_for_sensitivity" % j)
    m.setParam('OutputFlag', False)
    # m.params.IntFeasTol = 1e-7

    ## Project complete data,Project Tadeness,construction completion time
    CT = m.addVar(obj=0, vtype=GRB.CONTINUOUS, name="CT_%d" % j)

    ## Activity start time
    ST = {}
    project_activities = project_activity[project]
    for row in project_activities.nodes():
        ST[row] = m.addVar(obj=0, vtype=GRB.CONTINUOUS, name="ST_%d_%s" % (j, row))

    ## Review sequence z_ij
    ## move to annealing objective function

    m.update()

    ## Constrain 8: activity starting constrain
    # equation 20
    for a in project_activities.nodes():
        for r in project_activities.node[a]['resources']:
            m.addConstr(ST[a], GRB.GREATER_EQUAL, AT[j, r],
                        name="constraint_8_project_%d_activity_%s_resource_%s" % (j, a, r))

    ## Constrain 9 activity sequence constrain
    # equation 21
    for row1, row2 in project_activities.edges():
        m.addConstr(ST[row1] + project_activities.node[row1]['duration'], GRB.LESS_EQUAL,
                    ST[row2], name="constraint_9_project_%d_activity_%s_activity_%s" % (j, row1, row2))

    ## Constrain 10,11
    for row1 in project_activities.nodes():
        for row2 in project_activities.nodes():
            if row1 != row2 and len(list(
                    set(project_activities.node[row1]['rk_resources']).intersection(
                        project_activities.node[row2]['rk_resources']))) > 0:
                # equation 22
                m.addConstr(ST[row1] + project_activities.node[row1]['duration'] - M * (
                    1 - _get_y_for_activities(y_, row1, row2)), GRB.LESS_EQUAL, ST[row2],
                            name="constraint_10_project_%d_activity_%s_activity_%s" % (j, row1, row2))
                # equation 23
                m.addConstr(
                    ST[row2] + project_activities.node[row2]['duration'] - M * _get_y_for_activities(y_, row1, row2),
                    GRB.LESS_EQUAL, ST[row1],
                    name="constraint_11_project_%d_activity_%s_activity_%s" % (j, row1, row2))
                # m.addConstr(y[j,row1,row2]+y[j,row2,row1],GRB.LESS_EQUAL,1)

    ## Constrain 12
    # equation 24
    for row in project_activities.nodes():
        m.addConstr(CT, GRB.GREATER_EQUAL, ST[row] + project_activities.node[row]['duration'],
                    name="constraint_12_project_%d_activity_%s" % (j, row))

    m.update()

    # Set optimization objective - minimize completion time
    expr = LinExpr()
    expr.add(CT)
    m.setObjective(expr, GRB.MINIMIZE)
    m.update()
    ##########################################
    # m.params.presolve = 1
    m.update()
    m.setParam(GRB.Param.Method, 0)
    m.update()
    # Solve
    # m.params.presolve=0

    m.optimize()

    _skja = {}
    for c in m.getConstrs():
        if c.ConstrName.startswith('constraint_8_project'):
            splits = c.ConstrName.split('_')
            r = splits[7]
            if r not in _skja:
                _skja[r] = []
            _skja[r].append(c.Pi)

            # if c.Pi != 0:
            #     logging.debug('project %d binding resource:%s Pi:%.4g' % (j, splits[-1], c.Pi))
            # else:
            #     logging.debug('project %d not binding resource:%s Pi:%.4g' % (j, splits[-1], c.Pi))

            _pa_dataset.loc[_pa_dataset.shape[0]] = [_round, j, r, splits[5], c.Pi]
    _skj = {}
    for r in _skja:
        _skj[j, r] = max(_skja[r])
        _pa_max_dataset.loc[_pa_max_dataset.shape[0]] = [_round, j, r, max(_skja[r])]
    return _skj

# 需要导入模块: from gurobipy import Model [as 别名]
# 或者: from gurobipy.Model import getConstrs [as 别名]
def _sensitivity_analysis_for_tardiness(z, CT, D):
    m = Model("model_for_sensitivity_analysis_for_tardiness")
    m.setParam('OutputFlag', False)
    # m.params.IntFeasTol = 1e-7
    DT = {}
    TD = {}
    for j in range(D.project_n):
        ## Project Tadeness,construction completion time
        DT[j] = m.addVar(obj=0, vtype=GRB.CONTINUOUS, name="DT_%d" % j)
        TD[j] = m.addVar(obj=0, vtype=GRB.CONTINUOUS, name="TD_%d" % j)

    DT[-1] = m.addVar(obj=0, vtype=GRB.CONTINUOUS, name="DT_-1")
    m.update();

    #### Add Constraint ####
    ## Constrain 2: project complete data>due data ##
    # equation 17
    for j in range(D.project_n):
        m.addConstr(DT[j] - TD[j], GRB.LESS_EQUAL, D.DD[j], name="constraint_2_project_%d" % j)

    ## Constraint 13
    # equation 12
    for j in range(D.project_n):
        m.addConstr(DT[j], GRB.GREATER_EQUAL, CT[j] + D.review_duration[j],
                    name="constraint_13_project_%d" % j)

    ## Constraint 14
    # equation 13
    for i in range(-1, D.project_n):
        for j in range(D.project_n):
            if i != j:
                m.addConstr(DT[j], GRB.GREATER_EQUAL, DT[i] - D.M * (1 - z[i, j]) + D.review_duration[j],
                            name="constraint_14_project_%d_project_%d" % (i, j))

    m.update()

    # Set optimization objective - minimize sum of
    expr = LinExpr()
    for j in range(D.project_n):
        expr.add(D.w[j] * TD[j])
    m.setObjective(expr, GRB.MINIMIZE)
    m.update()

    # m.params.presolve = 1
    m.update()

    m.optimize()

    # _logger.info("mm binding info:")
    sj = {}
    for c in m.getConstrs():
        if c.ConstrName.startswith('constraint_13'):
            j = int(c.ConstrName.split('_')[-1])
            # if c.Pi != 0:
            # sj[j] = 1
            # _logger.info('%s binding Pi:%.4g' % (c.ConstrName, c.Pi))
            # pass
            # else:
            # sj[j] = 0
            # _logger.info('%s not binding Pi:%.4g' % (c.ConstrName, c.Pi))
            # pass
            sj[j] = c.Pi
    return sj

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

#.........这里部分代码省略.........
    ## move to annealing objective function

    # y
    y = {}
    for activity_i in project_activities.nodes():
        for activity_j in project_activities.nodes():
            # print(project_activities.node[activity_i])
            # print(dir(project_activities.node[activity_i]))
            if activity_i != activity_j and len(list(
                    set(project_activities.node[activity_i]['rk_resources']).intersection(
                        project_activities.node[activity_j]['rk_resources']))) > 0:
                y[activity_i, activity_j] = m.addVar(obj=0, vtype=GRB.BINARY,
                                                     name="(y%d,%s,%s)" % (j, activity_i, activity_j))
    m.update()

    #### Create constrains ####
    ## Constrain 2: project complete data>due data
    ## move to annealing objective function

    ## Constrain 3: supplier capacity limit
    ## move to annealing neighbor & random generator

    ## Constrain 4,6: project demand require; each project receive from one supplier for each resource
    ## move to annealing neighbor & random generator

    ## constrain 5: shipping constrain
    ## move to annealing neighbor & random generator

    ## Constrain 7:budget limit
    ## move to annealing constraint valid

    ## Constrain 8: activity starting constrain
    for a in project_activities.nodes():
        for r in project_activities.node[a]['resources']:
            m.addConstr(AT[j, r], GRB.LESS_EQUAL, ST[a],
                        name="constraint_8_project_%d_activity_%s_resource_%s" % (j, a, r))

    ## Constrain 9 activity sequence constrain
    for row1, row2 in project_activities.edges():
        m.addConstr(ST[row1] + project_activities.node[row1]['duration'], GRB.LESS_EQUAL,
                    ST[row2], name="constraint_9_project_%d_activity_%s_activity_%s" % (j, row1, row2))

    ## Constrain 10,11
    for row1 in project_activities.nodes():
        for row2 in project_activities.nodes():
            if row1 != row2 and len(list(
                    set(project_activities.node[row1]['rk_resources']).intersection(
                        project_activities.node[row2]['rk_resources']))) > 0:
                m.addConstr(ST[row1] + project_activities.node[row1]['duration'] - M * (
                    1 - y[row1, row2]), GRB.LESS_EQUAL, ST[row2],
                            name="constraint_10_project_%d_activity_%s_activity_%s" % (j, row1, row2))
                m.addConstr(
                    ST[row2] + project_activities.node[row2]['duration'] - M * (y[row1, row2]),
                    GRB.LESS_EQUAL, ST[row1],
                    name="constraint_11_project_%d_activity_%s_activity_%s" % (j, row1, row2))
                # m.addConstr(y[j,row1,row2]+y[j,row2,row1],GRB.LESS_EQUAL,1)

    ## Constrain 12
    for row in project_activities.nodes():
        m.addConstr(CT, GRB.GREATER_EQUAL, ST[row] + project_activities.node[row]['duration'],
                    name="constraint_12_project_%d_activity_%s" % (j, row))

    ## Constrain 13
    ## move to anealing objective function

    ## Constrain 14
    ## move to anealing objective function

    ## Constrain 15
    ## move to anealing objective function

    ## Constrain 16
    ## move to anealing objective function

    ## Constrain 17
    ## move to anealing objective function

    m.update()

    # Set optimization objective - minimize completion time
    expr = LinExpr()
    expr.add(CT)
    m.setObjective(expr, GRB.MINIMIZE)
    m.update()
    ##########################################
    m.params.presolve = 1
    m.update()
    # Solve
    # m.params.presolve=0

    m.optimize()
    # m.write(join(output_dir, "heuristic_%d.lp" % j))
    # m.write(join(output_dir, "heuristic_%d.sol" % j))
    # logging.info("project %d with optimalVal %r" % (j, m.objVal))
    print(m.status == GRB.OPTIMAL)
    for c in m.getConstrs():
        if c.ConstrName.startswith('constraint_8_project'):
            c.getAttr(GRB.Attr.SARHSLow)
            logging.info('%s shadow price (%f,%f)' % (c.ConstrName, c.SARHSLow, c.SARHSUp))
    return m.objVal