Python调用Gurobi求解简单线性规划问题

前言

Gurobi是一款功能强大的商用求解器,支持Python、C、C++、Java等多种语言调用,相比于Cplex,Gurobi封装更高,更加方便,但是对于初学者而言会更难一些。Gurobi与Cplex存在兼容,Gurobi生成的mps或者lp文件可以在Cplex中运行。另外,在校学生和教师可以免费使用Gurobi的学术版,求解规模和求解速度不受限制。本文将简单记录python调用Gurobi配置及其求解的过程。

安装

用户可以通过gurobi官网下载并安装软件:gurobi官网,下载软件后学生可通过学生认证获取免费证书,获取证书后在电脑上验证即可。获取license后,在终端运行:grbgetkey xxxxxxxx,如下图所示。

运行完成后会得到gurobi.lic文件,可以放在指定目录下。之后在anaconda对应的环境中安装gurobi包,安装语句为:

conda install -c gurobi gurobi

安装完成后,每次运行程序需要激活许可证,即上面提到的gurobi.lic文件,输入以下语句:

export GRB_LICENSE_FILE=/Users/liuxinwei/gurobi.lic

许可证验证完成后即可调用gurobi包进行数学规划等操作。

简单示例(一)

'''
线性规划模型
max x+y+2z
s.t.
    x+2y+3z<=4
    x+y>=1
    x,y,z=0 or 1
'''

from gurobipy import *

try:
    #模型
    model=Model('mip')

    #变量
    x=model.addVar(vtype=GRB.BINARY,name='x')
    y=model.addVar(vtype=GRB.BINARY,name='y')
    z=model.addVar(vtype=GRB.BINARY,name='z')

    #目标函数
    model.setObjective(x+y+2*z,GRB.MAXIMIZE)

    #约束
    model.addConstr(x+2*y+3*z<=4,name='c1')
    model.addConstr(x+y>=1,name='c2')

    #求解
    model.setParam('outPutFlag',0)#不输出求解日志
    model.optimize()

    #输出
    print('obj=',model.objVal)
    for v in model.getVars():
        print(v.varName,':',v.x)
        
except GurobiError as e:
    print('Error code '+str(e.errno)+':'+str(e))

except AttributeError:
    print('Encountered an attribute error')

简单示例(二)

'''
营养配方模型
目标:以最低的价格满足人体所需营养
约束:人体所需营养必须有上限和下限,每种食物有营养成分和价格的属性
'''

from gurobipy import *

#数据
categories,minNutrition,maxNutrition=multidict({
    'calories':[1800,2200],
    'protein':[91,GRB.INFINITY],
    'fat':[0,65],
    'sodium':[0,1779]
})

foods,cost=multidict({
    'hamburger':2.49,
    'chicken':2.89,
    'hot dog':1.50,
    'fries':1.89,
    'macaroni':2.09,
    'pizza':1.99,
    'salad':2.49,
    'milk':0.89,
    'ice cream':1.59
})

nutritionValues = {
    ('hamburger', 'calories'): 410,
    ('hamburger', 'protein'):  24,
    ('hamburger', 'fat'):      26,
    ('hamburger', 'sodium'):   730,
    ('chicken',   'calories'): 420,
    ('chicken',   'protein'):  32,
    ('chicken',   'fat'):      10,
    ('chicken',   'sodium'):   1190,
    ('hot dog',   'calories'): 560,
    ('hot dog',   'protein'):  20,
    ('hot dog',   'fat'):      32,
    ('hot dog',   'sodium'):   1800,
    ('fries',     'calories'): 380,
    ('fries',     'protein'):  4,
    ('fries',     'fat'):      19,
    ('fries',     'sodium'):   270,
    ('macaroni',  'calories'): 320,
    ('macaroni',  'protein'):  12,
    ('macaroni',  'fat'):      10,
    ('macaroni',  'sodium'):   930,
    ('pizza',     'calories'): 320,
    ('pizza',     'protein'):  15,
    ('pizza',     'fat'):      12,
    ('pizza',     'sodium'):   820,
    ('salad',     'calories'): 320,
    ('salad',     'protein'):  31,
    ('salad',     'fat'):      12,
    ('salad',     'sodium'):   1230,
    ('milk',      'calories'): 100,
    ('milk',      'protein'):  8,
    ('milk',      'fat'):      2.5,
    ('milk',      'sodium'):   125,
    ('ice cream', 'calories'): 330,
    ('ice cream', 'protein'):  8,
    ('ice cream', 'fat'):      10,
    ('ice cream', 'sodium'):   180}

#模型
model=Model('diet')

#变量
buy=model.addVars(foods,name='buy')

#目标函数
model.setObjective(buy.prod(cost),GRB.MINIMIZE)
#model.setObjective(sum(buy[f]*cost[f] for f in foods),GRB.MINIMIZE)

#约束(expr==[lb,ub]-->lb<=expr<=ub)
model.addConstrs((quicksum(nutritionValues[f,c]*buy[f] for f in foods)\
                 ==[minNutrition[c],maxNutrition[c]] for c in categories),'con')
#求解
model.setParam('OutputFlag', 0)
model.optimize()

#输出
if model.status==GRB.Status.OPTIMAL:
    print('obj=',model.objVal)
    buyx=model.getAttr('x',buy)#取buy的x属性
    for f in foods:
        if buyx[f]>0.0001:
            print(f,':',buyx[f])

简单示例(三)

'''
指派模型
目标:成本最小化,时间最小化
'''

from gurobipy import *
import random

#数据
N=3
random.seed(1)
T={(i,j):random.randint(0,100) for i in range(N) for j in range(N)}
C={(i,j):random.randint(0,100) for i in range(N) for j in range(N)}

#模型
model=Model('assignment')

#添加变量
x=model.addVars([(i,j) for i in range(N) for j in range(N)],vtype='B',name='x')

#目标函数
model.setObjectiveN(x.prod(T),index=0,weight=0.1,name='Time Objective')
model.setObjectiveN(x.prod(C),index=1,weight=0.5,name='Cost Objective')

#约束
model.addConstrs(((quicksum(x[i,j] for i in range(N))==1) for j in range(N)),name='Job Constraingt')
model.addConstrs(((quicksum(x[i,j] for j in range(N))==1) for i in range(N)),name='Worker Constraingt')

#求解
model.Params.OutputFlag=0
model.optimize()

#输出
for i in range(model.NumObj):
    model.setParam('ObjNumber',i)
    print('Obj{}:'.format(i),model.objVal)

for i in range(N):
    for j in range(N):
        if x[i,j].x>0:
            print('job {}-->worker {}'.format(j+1,i+1))
            
model.write('assignment.lp')

Callback函数调用

'''
callback函数调用
'''

from gurobipy import *
import random

#callback函数
def RINScallback(model,where):
    if where==GRB.Callback.MIPNODE:
        #节点数量%100==0时,调用callback函数
        if model.cbGET(GRB.Callback.MINPNODE_NODCNT)%100==0 and\
           model.cbGET(GRB.Callback.MIPNODE.STATUS)==GRB>OPTIMAL:
            submodel=model.copy()
            suby=submodel.getVars()
            
            #获得节点松弛解
            yrelaxation=model.cbGetNodeRel(model._y)
            
            #固定变量取值
            for i in range(model._N):
                if abs(yrelaxation[i])<0.01:
                    suby[i].ub=0
                elif abs(yrelaxation-1)<=0.01:
                    suby[i].lb=1
            submodel.setParam(GRB.Param.TimeLimit,30)
            submodel.optimize()
            if submodel.objval>model.cbGet(GRB.Callback.MIPNODE_OBJBST):
                #将解传递给原模型
                for i in range(model._N):
                    if abs(suby[i].x)<=0.001:
                        model.cbSetSolution(model._y[i],0.0)
                    elif abs(suby[i].x-1)<=0.001:
                        model.cbSetSolution(model._y[i],1.0)

#数据
random.seed(1)
N=10
C={(i,j):random.randint(0,100) for i in range(N) for j in range(N)}

#模型
model=Model('MaxCut')

#变量
y=model.addVars(N,vtype='B',name='y')

#目标函数
obj=QuadExpr()
for i in range(N):
    for j in range(N):
        obj+=C[i,j]*(2*y[i]-1)*(2*y[j]-1)
model.setObjective(obj,GRB.MAXIMIZE)

#设置求解时间
model.Params.TimeLimit=600

#外部变量
model._y=y
model._N=N

#求解
model.optimize(RINScallback)

广义约束及其线性化

#模型
model=Model('gurobi')

#变量
x=model.addVar(vtype='C',name='x')
y=model.addVar(vtype='C',name='y')
z=model.addVar(vtype='C',name='z')

#目标
model.setObjective(x+y+z,GRB.MAXIMIZE)

#约束
model.addConstr(x+y<=10)
model.addConstr(z==min_(x,y,6))

#求解
model.Params.OutputFlag=0
model.optimize()

#输出
print('obj=',model.objVal)
for v in model.getVars():
    print(v.varName+'='+str(v.x))

另外,对于取绝对值的情况,例如y=|x|,可写成model.addConstr(y==abs_(x))。


如上约束也可写成:
min z
x<=My
y={0,1}
z>=cx+k+(y-1)M

对于分段式目标函数,我们也可通过变量的转化将其进行线性化处理:

对于逻辑或的约束,即两个约束至少有一个成立的问题,可按照如下方式进行线性化:

总结

后续将继续跟进Gurobi的学习,做一些复杂问题的建模求解。

物联沃分享整理
物联沃-IOTWORD物联网 » Python调用Gurobi求解简单线性规划问题

发表评论