In [1]:
#線形最適化問題
from gurobipy import *
In [21]:
#モデルを作成
model = Model("k")
In [22]:
#modelに変数を追加 addVar
x1 = model.addVar(ub =5, lb=0 ,name="x1")
x2 = model.addVar(ub =6 ,lb=0, name="x2")
x3 = model.addVar(ub = 4, lb=0,name="x3")
model.update()
In [23]:
#制約の追加 addConstr
model.addConstr(x1 + x2 + x3 == 7)
Out[23]:
In [24]:
#目的関数の設定 setObjective
model.setObjective(121- x1 - 3*x2 - 4*x3)
In [25]:
#目的関数の方向を設定 ModelSense = -1で最大化 1で最小化 規定値は1
model.ModelSense =1
In [26]:
#解の出力 ObjVal
from __future__ import division
model.optimize()
print "Opt. Value= " ,model.ObjVal
In [27]:
#会の出力その2 getVarsは変数オブジェクトのリストを返す VarNameとX
for v in model.getVars():
print v.VarName, v.X
In [24]:
#これでも最適解を出力可能
print "(x1, x2 ,x3)=" , x1.X, x2.X, x3.X
In [25]:
#整数最適化問題 ほとんど上と同じ vtypeをIにするだけ
model = Model("turukame")
x = model.addVar(vtype = "I")
y = model.addVar(vtype = "I")
z = model.addVar(vtype = "I")
model.update()
model.addConstr(x + y+ z == 32)
model.addConstr(2*x + 4*y +8*z == 80)
model.setObjective(y + z )
model.optimize()
print "Opt. Value =" ,model.ObjVal
print "(x,y,z)=" , x.X ,y.X ,z.X
In [72]:
#輸送問題
model = Model()
In [73]:
#顧客リストと需用量辞書、工場リストと供給量リストの同時出力 multidict
I,d = multidict({1:80, 2:270 ,3:250 ,4:160 ,5:180})
J,M = multidict({1:500 ,2:500 ,3:500})
In [35]:
I
Out[35]:
In [37]:
J
Out[37]:
In [38]:
d
Out[38]:
In [39]:
M
Out[39]:
In [74]:
#輸送費用の設定 タプルで顧客と工場の関係を表示 タプルをキーとする辞書型
c = {(1,1):4, (1,2):6, (1,3):9, (2,1):5, (2,2):4, (2,3):7, (3,1):6, (3,2):3, (3,3):4, (4,1):8, (4,2):5, (4,3):3, (5,1):10, (5,2):8, (5,3):4}
In [40]:
c
Out[40]:
In [75]:
#変数の作成
x = {}
for i in I:
for j in J:
x[i,j] = model.addVar(vtype = "C", name = "x(%s,%s)" %(i,j))
model.update()
In [21]:
x
Out[21]:
In [76]:
#制約の追加 オプションで制約にも名前をつけていることに注意 ここではそれぞれDeman(i)という制約にしている。
for i in I:
model.addConstr(quicksum(x[i,j] for j in J) == d[i], name = "Demand(%s)" % i)
In [77]:
for j in J:
model.addConstr(quicksum(x[i,j] for i in I ) <= M[j] ,name = "Capacity(%s)"% j )
In [29]:
#次の二つはどちらでも良い
c[(1,2)]
Out[29]:
In [30]:
c[1,2]
Out[30]:
In [78]:
#目的関数の設定
model.setObjective(quicksum(c[i,j]*x[i,j] for (i,j) in x) ,GRB.MINIMIZE)
In [79]:
model.optimize()
In [81]:
#結果の表示
print "optimal value:" , model.ObjVal
EPS = 1.e-6 #0以下の結果を表示しないために
for (i,j) in x:
if x[i,j].X > EPS:
print "sending quantity %10s from %3s to customer %3s" %(x[i,j].X,j,i) #Xで値だけ呼ばないと、編すの持つ全ての要素が呼び出される。
In [83]:
#双対問題
print "Const. Name : Slack ,Dual"
for c in model.getConstrs():
print "%s : %s, %s" %(c.ConstrName , c.Slack, c.pi)
In [84]:
#混合問題
#各原料が含む成分の割合
a = {(1,1):.25 ,(1,2):.15 ,(1,3):.3, (2,1):.3 ,(2,2):.3 ,(2,3):.1 ,(3,1):.15 ,(3,2):.65 ,(3,3):.05 ,(4,1):.1 ,(4,2):.05 ,(4,3):.85}
In [88]:
#各種データ
I,p = multidict({1:5, 2:6 ,3:8 ,4:20})
K,LB,UB = multidict({1:[.1,.2], 2:[.0,.35], 3:[.45,1.0]})
In [90]:
model = Model("product_mix")
x ={}
for i in I:
x[i] = model.addVar()
model.update()
model.addConstr(quicksum(x[i] for i in I) == 1)
for k in K:
model.addConstr(quicksum(a[i,k] * x[i] for i in I ) <= UB[k])
model.addConstr(quicksum(a[i,k] * x[i] for i in I ) >= LB[k])
model.setObjective(quicksum(p[i] * x[i] for i in I))
model.optimize()
for i in I:
print i, x[i].X
In [2]:
#分数最適化 目的関数を制約に取り込んで二分探索法
LB, UB, EPS = 0.0 ,1.0 ,0.01
while 1:
theta = (UB+LB)/2
model = Model("tako")
x = model.addVar(vtype = "I")
y = model.addVar(vtype = "I")
z = model.addVar(vtype = "I")
model.update()
model.addConstr(x+y+z ==32)
model.addConstr(2*x +4*y +8*z == 80)
model.addConstr((2*theta-1)*x +(4*theta-1)*y >=0)
model.setObjective(x+y+z)
model.optimize()
if model.Status == GRB.OPTIMAL:
UB = theta
if UB - LB <= EPS:
break
else:
LB = theta
print "(x,y,z)=" , x.X ,y.X ,z.X
In [3]:
#多制約0-1ナップザック問題
#まず問題の定式化と必要なデータ型の準備をする
def example():
J,v = multidict({1:16,2:19,3:23,4:28})
a = {(1,1):2 ,(1,2):3 ,(1,3):4 ,(1,4):5 ,(2,1):3000 ,(2,2):3500 ,(2,3):5100 ,(2,4):7200}
I,b = multidict({1:7, 2:10000})
return I ,J ,v ,a,b
In [4]:
def mkp(I ,J ,v ,a ,b):
model = Model("nap")
x ={}
for j in J:
x[j] = model.addVar(vtype ="B" ,name ="x(%d)"%j) #%dは符号付整数 ブール型
model.update()
for i in I:
model.addConstr(quicksum(a[i,j]*x[j] for j in J) <= b[i])
model.setObjective(quicksum(v[j] * x[j] for j in J), GRB.MAXIMIZE)
model.update()
return model
In [5]:
#最適化本番
if __name__ == "__main__": #最適化本番に入る際のおまじない
I ,J ,v ,a ,b = example()
model = mkp(I ,J ,v ,a ,b)
In [6]:
#定式化のデバック（ただし、デバック画面を見れない）
model.update()
model.write("mkp.lp")
In [8]:
#最適化
model.optimize()
In [9]:
print "Optimum .Value=" ,model.ObjVal
EPS = 1.e-6
for v in model.getVars():
if v.X > EPS:
print v.VarName ,v.X
In [21]:
#栄養問題
model = Model("nutrients")
F =["c", "b", "d"]
x ={}
for j in F:
x[j] = model.addVar(vtype ="I" ,name="x(%s)" %j) #%sは文字列への変換だから納得
model.update()
In [24]:
#規約不整合部分系 computeIIS
#逸脱最適化 feasRelaxS
F,c,d = multidict({"CQPounder": [ 360, {"Cal":556, "Carbo":39, "Protein":30,"VitA":147,"VitC": 10, "Calc":221, "Iron":2.4}],
"Big Mac" : [ 320, {"Cal":556, "Carbo":46, "Protein":26,
"VitA":97, "VitC": 9, "Calc":142, "Iron":2.4}],
"FFilet" : [ 270, {"Cal":356, "Carbo":42, "Protein":14,
"VitA":28, "VitC": 1, "Calc": 76, "Iron":0.7}],
"Chicken" : [ 290, {"Cal":431, "Carbo":45, "Protein":20,
"VitA": 9, "VitC": 2, "Calc": 37, "Iron":0.9}],
"Fries" : [ 190, {"Cal":249, "Carbo":30, "Protein": 3,
"VitA": 0, "VitC": 5, "Calc": 7, "Iron":0.6}],
"Milk" : [ 170, {"Cal":138, "Carbo":10, "Protein": 7,
"VitA":80, "VitC": 2, "Calc":227, "Iron": 0}],
"VegJuice" : [ 100, {"Cal": 69, "Carbo":17, "Protein": 1,
"VitA":750,"VitC": 2, "Calc":18, "Iron": 0}],})
N,a,b = multidict({
"Cal" : [ 2000, 3000],
"Carbo" : [ 300, 375 ],
"Protein" : [ 50, 60 ],
"VitA" : [ 500, 750 ],
"VitC" : [ 85, 100 ],
"Calc" : [ 660, 900 ],
"Iron" : [ 6.0, 7.5 ],
})
In [26]:
d
Out[26]:
In [6]:
#なぜかinfeasible
model2 =Model("diet")
x ={}
for j in F:
x[j] = model2.addVar(vtype="I" ,name = "x(%s)" %j)
model2.update()
for i in N:
model2.addConstr(quicksum(d[j][i]*x[j] for j in F) >= a[i] ,name ="NutrLB(%s)"%i)
model2.addConstr(quicksum(d[j][i]*x[j] for j in F) <= b[i] ,name = "NutrUB(%s)"%i)
model2.setObjective(quicksum(c[j] *x[j] for j in F))
model2.optimize()
In [ ]: