カーネルやら何やらを理解するために、非線形計画法を勉強中。
で、その前に線形計画法で実際に手を動かしたくなったので無償で使えるPuLPを使ってみた。
今回例題として解いたのは「これならわかる最適化数学」の第6章、線形計画法。
似たような例題がいくつか並ぶので、ピックアップしてPuLPでモデル化してみた。
【例題6.2】(P161)
2種類の容器A,Bを作るのに, 機械M1,M2を使う. 容器Aを1個作るのにM1を2分, 機械M2を4分使う必要がある.一方,容器Bを1個作るのに、機械M1を8分,機械M2を4分使う必要がある. 容器A,Bを作る利益は一つあたりそれぞれ29円,45円である. 利益を最大にするにはどのように計画すれば良いか.
利益を目的関数f、容器の個数をそれぞれx,yを問いてモデル化すると式は下記のようになる。
制約条件:
2x + 8y <= 60
4x + 4y <= 60
x >= 0, y >= 0
目的関数：
f = 29x + 45y → max
で、これらの式をそのまんまPuLPに突っ込んで解かせるコードが下記。
今回は容器の個数を求めるので、解は整数であり、変数はIntegerで指定する。
#Creat variables
x = pulp.LpVariable("x", cat = "Integer")
y = pulp.LpVariable("y", cat = "Integer")
#Create object function
problem = pulp.LpProblem("Container", pulp.LpMaximize)
problem += 29*x + 45*y #(6.6)
#constrained condition
problem += 2*x + 8*y <= 60 #(6.4)
problem += 4*x + 4*y <= 60 #(6.4)
problem += x >= 0 #(6.5)
problem += y >= 0 #(6.5)
status = problem.solve()
print "Status", pulp.LpStatus[status]
print problem
print "Result"
print "x", x.value()
print "y", y.value()
結果は、
Status Optimal
Container:
MAXIMIZE
29*x + 45*y + 0
SUBJECT TO
_C1: 2 x + 8 y <= 60
_C2: 4 x + 4 y <= 60
_C3: x >= 0
_C4: y >= 0
VARIABLES
x free Integer
y free Integer
Result
x 10.0
y 5.0
【例題6.3】(P161)
2種類の金属M1,M2を用いて2種類の合金A,Bを作る. 合金A,Bはそれぞれ1トン当たり30,000円、25,000円の利益がある。合金Aは金属M1,M2を1:1の割合で混合し、合金Bは金属M1,M2を1:3の割合で混合する。金属M1,M2はそれぞれ1日当たり10トン,15トンの供給が可能である. 利益を最大にするのに合金A,Bを作るにはどのように計画すればよいか.
モデル化すると、
制約条件:
0.5x + 0.25y <= 10
0.5x + 0.75y <= 15
x >= 0, y >= 0
目的関数：
f = 30x + 25y → max
今回はコード化にあたり、解は連続値であるので、変数はContinuousで指定する。
#Creat variables
x = pulp.LpVariable("x", cat = "Continuous")
y = pulp.LpVariable("y", cat = "Continuous")
#Create object function
problem = pulp.LpProblem("Alloy", pulp.LpMaximize)
problem += 30*x + 25*y #(6.9)
#constrained condition
problem += 0.5*x + 0.25*y <= 10 #(6.7)
problem += 0.5*x + 0.75*y <= 15 #(6.7)
problem += x >= 0 #(6.8)
problem += y >= 0 #(6.8)
status = problem.solve()
print "Status", pulp.LpStatus[status]
print problem
print "Result"
print "x", x.value()
print "y", y.value()
結果は、
Status Optimal
Alloy:
MAXIMIZE
30*x + 25*y + 0
SUBJECT TO
_C1: 0.5 x + 0.25 y <= 10
_C2: 0.5 x + 0.75 y <= 15
_C3: x >= 0
_C4: y >= 0
VARIABLES
x free Continuous
y free Continuous
Result
x 15.0
y 10.0
【例題6.7】(P164)
今回の例のように目的関数が発散して、最適解が存在しない場合は下記のように、"Status Unbounded"と表示される。
制約条件:
-x - y <= -1
-2x + y <= 1
x - 2y <= 1
x >= 0, y >= 0
目的関数：
f = x + y → max
#Creat variables
a = pulp.LpVariable("a", cat = "Continuous")
b = pulp.LpVariable("b", cat = "Continuous")
#Create object function
problem = pulp.LpProblem("Test", pulp.LpMaximize)
problem += a + b #(6.16)
#constrained condition
problem += -a - b <= -2 #(6.14)
problem += -2*a + b <= 15 #(6.14)
problem += a - 2*b <= 15 #(6.14)
problem += a >= 0 #(6.15)
problem += b >= 0 #(6.15)
status = problem.solve()
print "Status", pulp.LpStatus[status]
print problem
Status Unbounded