（R Advent Calendar @ Qiita 16日目）
こんにちは、16日目担当のサーモンさんです。
突然なんですが、僕Linear Programming好きなんですよ。
というか多分Operations ResearchとかSupply Chain Managementが好きなんですよ。
なんなんすかね、ビジネスのフローをモデル化する部分が好きなんでしょう。
卒論でもサーモン養殖の線形計画法とかやろうか悩んでたくらいですよ。（収穫船と餌の輸送と生産の組合せ計画とかがあるんすよ）
で、大学院にいた時はAMPLっていう言語でCPLEXとかNEOSとかっていったソルバーを使って問題を解いていたのですが、ソルバーはお高い訳です。
まぁ早い話がRでやりたい訳です。
で、探してみる訳ですよ。
で、理論とかのお話が多く、あんまり例題が無い訳です。
という事で今回は線形計画法の応用についてを書きたいと思います。
1. 生産計画問題(Production Planning)
ある小さな会社は3つのタイプのプロダクトを生産しています。
Model Sがスタンダードなタイプで、Model LとPはそれに多少の変更をしたタイプで、これらの商品は(punching, forming, assembly)の3つの生産過程を経て完成されます。
それぞれの生産過程では配置されている人員が固定されている為に、労働量がそれぞれの過程に置いて限界があります。また、3つのプロダクトはそれぞれ生産に掛かる労力が異なっています。下の表がこれらの条件をまとめた物となっています。(AssemblyにおいてはSとLPが別のプロセスになっています。)
さて、3つのプロダクト(S,L,P)の利益がそれぞれ(12o,190,180)だった時、どのような組み合わせで生産を行えば利益が最大化出来るのか？
というのがここで考えたい事です。
xをプロダクトの生産量として、モデルを組んでみると以下の様になります。
1行目は目的関数で、これを最大化します。中身としては(利益＊生産量)をそれぞれのプロダクトで計算して合計すると言う物です。
2行目以降が変数の制約条件となっています。
2,3,4,5行目は生産過程の労働量の制約です。中身としては、(1生産量当たりに掛かる労力*生産量)をプロダクトごとに出して合計し、それが用意されている労働量(Capacity)を超えない様な状態になる様に指定しています。
6行目は生産量が負にならない様に制約をしています。
さて、モデルを組んだので後はRの出番です。
library(lpSolveAPI)
#define model
lpmodel<-make.lp(4,3)
#define variable and parameter
set.column(lpmodel,1,c(120,0.3,0.25,0.3), indices=c(0,1,2,3))
set.column(lpmodel,2,c(190,0.3,0.5,0.6), indices=c(0,1,2,4))
set.column(lpmodel,3,c(180,0.3,0.45,0.5), indices=c(0,1,2,4))
#define constraints
set.constr.value(lpmodel, rhs=800, constraints=1)
set.constr.value(lpmodel, rhs=800, constraints=2)
set.constr.value(lpmodel, rhs=570, constraints=3)
set.constr.value(lpmodel, rhs=840, constraints=4)
#max or min
lp.control(lpmodel,sense=’max’)
#write model file
write.lp(lpmodel,’model.lp’,type=’lp’)
#solve the model, if this return 0 an optimal solution is found
solve(lpmodel)
#this return the proposed solution
get.objective(lpmodel)
みての通り今回利用したパッケージはlpSolveAPIです。
理由は特に無いですw
パッケージのリファレンスを見た時に必要そうな機能が全部あった事と、モデルが自分としては書きやすかった事ですかね。
後はソルバーに渡す前のモデルファイルを自分でちゃんと見れるのが結構便利かなと。
さて、やっている事の解説をして行きます。
まず最初に変数の数と制約条件の数を指定してモデルの枠組みを作ります。
今回はプロダクトが3つで、製造過程での制約が4つなので、以下の様に書きます。
lpmodel<-make.lp(4,3)
これでモデルの枠組みが作られました。
この時モデルは以下の様な構造になっています。
次に変数とパラメーターの設定を行います。
set.column(lpmodel,1,c(120,0.3,0.25,0.3), indices=c(0,1,2,3))
set.column(lpmodel,2,c(190,0.3,0.5,0.6), indices=c(0,1,2,4))
set.column(lpmodel,3,c(180,0.3,0.45,0.5), indices=c(0,1,2,4))
定義としてはこんな感じです。
set.column(モデル名, Columnの指定, 変数に付くパラメーター, indiciesの指定)
例えば一番上のを実行すると、変数1の目的関数と制約条件1,2,3に、(120*変数1),(0.3*変数1),(0.25*変数1),(0.3*変数1)を入れてくれます。
2番目を実行すると、変数2の目的関数と制約条件1,2,4に、(190*変数2),(0.3*変数2),(0.5*変数2),(0.6*変数2)を入れてくれます。
3番目を実行すると、変数3の目的関数と制約条件1,2,4に、(180*変数2),(0.3*変数2),(0.45*変数3),(0.5*変数3)を入れてくれます。
結果こんな感じになります。列の間は+で繋がれているるので、例えば目的関数であれば以下の式が出来ている事になります。
120*x1 + 190*x2 + 180*x3
次に制約条件の右辺を入れてゆきます。
set.constr.value(lpmodel, rhs=800, constraints=1)
set.constr.value(lpmodel, rhs=800, constraints=2)
set.constr.value(lpmodel, rhs=570, constraints=3)
set.constr.value(lpmodel, rhs=840, constraints=4)
set.constr.value(モデル名, rhs=右辺の値, constraints=制約条件の番号（indicies))
これでさっきの表の右辺に値が入ってゆきます。
デフォルトでは左辺<=右辺という設定がされているので、今回は特に左辺と右辺の関係を調整する必要はありません。
最後に目的が最小化なのか、最大化なのかを定義します。
lp.control(lpmodel,sense=’max’)
後はモデルを書いて、実行するだけです。
write.lp(lpmodel,’model.lp’,type=’lp’)
これで上で定義してきたモデルを書き込みます。（別段必要ない気もしますが。）
モデルファイルはテキストで開けるので、エラーが起きた際なんかに自分の書いてるモデルがどんなもんなのかが確認できて良いです。
今回作ったモデルは下の様な形で記述されています。
/* Objective function */
max: +120 C1 +190 C2 +180 C3;
/* Constraints */
+0.3 C1 +0.3 C2 +0.3 C3 <= 800;
+0.25 C1 +0.5 C2 +0.45 C3 <= 800;
R3: +0.3 C1 <= 570;
+0.6 C2 +0.5 C3 <= 840;
solve(lpmodel)
これでモデルの最適化を解かせます。実行結果が0ならば解が見つかったことを指します。
get.objective(lpmodel)
これで解いた階の目的関数の値が出てきます。今回で言えば利益の総量です。
> get.objective(lpmodel)
[1] 358000
この時の変数の値、つまり各商品の生産量はget.variables()で出てきます。
> get.variables(lpmodel)
[1] 1900.0000 0.0000 722.2222
つまり、プロダクトS = 1900, L = 0, P = 722.2222という事です。
2. 輸送計画問題(Transport and distribution problem)
次は輸送計画の問題です。
とある牛乳会社は4つの生産施設と5つの顧客を有しています。
それぞれの生産施設(F)は生産量に制限があり、5箇所の需要(K)もそれぞれ特定の量が決まっています。(単位：1000)
また各FとKの間には特定の輸送ルートが存在し、運搬量に応じて運送コストが掛かってしまいます。
これらの値は以下の表のとおりです。
ここで会社が行いたいことは、需要を満たしつつコストが最小限になるような生産と配送の組み合わせを見つける事です。
これをモデル化すると以下のようになります。（数列の表現にしています）
目的関数は運送コストの合計です。
一つ目の制約条件は、ある生産施設から5つの需要に対して送られる運送量が生産容量を超えないようにするものです。
二つ目の制約条件は、ある需要に対して各生産施設に送られてくる量の合計が需要量と合致しているように指定するものです。
これをRで解いてみると以下のような感じになります。
lpmodel<-make.lp(9,20)
set.column(lpmodel,1,c(2.8,1,1), indices=c(0,1,5))
set.column(lpmodel,2,c(2.55,1,1), indices=c(0,1,6))
set.column(lpmodel,3,c(3.25,1,1), indices=c(0,1,7))
set.column(lpmodel,4,c(4.3,1,1), indices=c(0,1,8))
set.column(lpmodel,5,c(4.35,1,1), indices=c(0,1,9))
set.column(lpmodel,6,c(4.3,1,1), indices=c(0,2,5))
set.column(lpmodel,7,c(3.15,1,1), indices=c(0,2,6))
set.column(lpmodel,8,c(2.55,1,1), indices=c(0,2,7))
set.column(lpmodel,9,c(3.3,1,1), indices=c(0,2,8))
set.column(lpmodel,10,c(3.5,1,1), indices=c(0,2,9))
set.column(lpmodel,11,c(3,1,1), indices=c(0,3,5))
set.column(lpmodel,12,c(3.3,1,1), indices=c(0,3,6))
set.column(lpmodel,13,c(2.9,1,1), indices=c(0,3,7))
set.column(lpmodel,14,c(4.3,1,1), indices=c(0,3,8))
set.column(lpmodel,15,c(3.4,1,1), indices=c(0,3,9))
set.column(lpmodel,16,c(5.2,1,1), indices=c(0,4,5))
set.column(lpmodel,17,c(4.45,1,1), indices=c(0,4,6))
set.column(lpmodel,18,c(3.5,1,1), indices=c(0,4,7))
set.column(lpmodel,19,c(3.75,1,1), indices=c(0,4,8))
set.column(lpmodel,20,c(2.45,1,1), indices=c(0,4,9))
set.constr.value(lpmodel, rhs=30000, constraints=1)
set.constr.value(lpmodel, rhs=40000, constraints=2)
set.constr.value(lpmodel, rhs=30000, constraints=3)
set.constr.value(lpmodel, rhs=40000, constraints=4)
set.constr.value(lpmodel, rhs=20000, constraints=5)
set.constr.value(lpmodel, rhs=30000, constraints=6)
set.constr.value(lpmodel, rhs=15000, constraints=7)
set.constr.value(lpmodel, rhs=25000, constraints=8)
set.constr.value(lpmodel, rhs=20000, constraints=9)
set.constr.type(lpmodel,rep(“=”,5),5:9)
lp.control(lpmodel,sense=’min’)
write.lp(lpmodel,’model.lp’,type=’lp’)
#solve the model, if this return 0 an optimal solution is found
solve(lpmodel)
#this return the proposed solution
get.objective(lpmodel)
生産施設が4つで需要が5つなので、合計２０個の意思決定の変数があります。
よってset.column()で２０個変数を作ります。
生産側で生産量の上限の制約があり、需要側で需要量と受け取る量を同じにする制約があるので、合計９個の制約条件があることになります。
よってset.constraint()で9個の制約を作ります。
そして、今回は需要と受け取りを=で結ぶ必要があるので、制約条件の5,6,7,8,9に関しては
set.constr.type(lpmodel, rep(“=”,5) , 5:9)
で制約条件の種類を変更しています。
実行結果は以下の通りです。
> get.objective(lpmodel)
[1] 306250
> get.variables(lpmodel)
[1] 0 30000 0 0 0 0 0 15000 25000 0 20000 0 0 0 0 0 0 0 0 20000
3. Transshipment problem
さて、実際の運送問題のケースではどこかに物資を中継させるというケースも存在します。
例えばこんな感じですね。先ほどの配送先に中継所D1,D2を加えてみます。
全ての生産施設はこの中継所D1,D2に以下の費用で牛乳を送ることが出来ます。（上の図では簡略化してますが。）
また、これら二つの中継所はすべての需要に以下の費用で牛乳を送ることが出来ます。(単位：1000)
これらを追加したモデルは以下のようになります。
生産と配送料の制約条件では、中継所への運送量も左辺に加わります。
二つ目の制約条件は、ある中継所に集まってきた牛乳がすべて送り出されるような制約を付けています。
３つ目の制約条件は、生産施設からの直接の運送と、中継所からの運送の合計が需要量とイコールになるようにしています。
これをRで解くとこんな感じになります。
lpmodel<-make.lp(11,38)
#x(i,j)
#i = 1
set.column(lpmodel,1,c(2.8,1,1), indices=c(0,1,5))
set.column(lpmodel,2,c(2.55,1,1), indices=c(0,1,6))
set.column(lpmodel,3,c(3.25,1,1), indices=c(0,1,7))
set.column(lpmodel,4,c(4.3,1,1), indices=c(0,1,8))
set.column(lpmodel,5,c(4.35,1,1), indices=c(0,1,9))
#i = 2
set.column(lpmodel,6,c(4.3,1,1), indices=c(0,2,5))
set.column(lpmodel,7,c(3.15,1,1), indices=c(0,2,6))
set.column(lpmodel,8,c(2.55,1,1), indices=c(0,2,7))
set.column(lpmodel,9,c(3.3,1,1), indices=c(0,2,8))
set.column(lpmodel,10,c(3.5,1,1), indices=c(0,2,9))
#i = 3
set.column(lpmodel,11,c(3,1,1), indices=c(0,3,5))
set.column(lpmodel,12,c(3.3,1,1), indices=c(0,3,6))
set.column(lpmodel,13,c(2.9,1,1), indices=c(0,3,7))
set.column(lpmodel,14,c(4.3,1,1), indices=c(0,3,8))
set.column(lpmodel,15,c(3.4,1,1), indices=c(0,3,9))
#i = 4
set.column(lpmodel,16,c(5.2,1,1), indices=c(0,4,5))
set.column(lpmodel,17,c(4.45,1,1), indices=c(0,4,6))
set.column(lpmodel,18,c(3.5,1,1), indices=c(0,4,7))
set.column(lpmodel,19,c(3.75,1,1), indices=c(0,4,8))
set.column(lpmodel,20,c(2.45,1,1), indices=c(0,4,9))
#y(i,k)
#k = 1
set.column(lpmodel,21,c(2.725,1,1), indices=c(0,1,10))
set.column(lpmodel,22,c(2.675,1,1), indices=c(0,2,10))
set.column(lpmodel,23,c(2.2,1,1), indices=c(0,3,10))
set.column(lpmodel,24,c(2.7,1,1), indices=c(0,4,10))
#k = 2
set.column(lpmodel,25,c(2.275,1,1), indices=c(0,1,11))
set.column(lpmodel,26,c(2.575,1,1), indices=c(0,2,11))
set.column(lpmodel,27,c(2.275,1,1), indices=c(0,3,11))
set.column(lpmodel,28,c(3.075,1,1), indices=c(0,4,11))
#w(k,j)
#k = 1
set.column(lpmodel,29,c(0.95,1,-1), indices=c(0,5,10))
set.column(lpmodel,30,c(0.75,1,-1), indices=c(0,6,10))
set.column(lpmodel,31,c(0.325,1,-1), indices=c(0,7,10))
set.column(lpmodel,32,c(0.95,1,-1), indices=c(0,8,10))
set.column(lpmodel,33,c(0.45,1,-1), indices=c(0,9,10))
#k = 2
set.column(lpmodel,34,c(0.7,1,-1), indices=c(0,5,11))
set.column(lpmodel,35,c(0.3,1,-1), indices=c(0,6,11))
set.column(lpmodel,36,c(0.3,1,-1), indices=c(0,7,11))
set.column(lpmodel,37,c(0.9,1,-1), indices=c(0,8,11))
set.column(lpmodel,38,c(0.875,1,-1), indices=c(0,9,11))
set.constr.value(lpmodel, rhs=30000, constraints=1)
set.constr.value(lpmodel, rhs=40000, constraints=2)
set.constr.value(lpmodel, rhs=30000, constraints=3)
set.constr.value(lpmodel, rhs=40000, constraints=4)
set.constr.value(lpmodel, rhs=20000, constraints=5)
set.constr.value(lpmodel, rhs=30000, constraints=6)
set.constr.value(lpmodel, rhs=15000, constraints=7)
set.constr.value(lpmodel, rhs=25000, constraints=8)
set.constr.value(lpmodel, rhs=20000, constraints=9)
set.constr.value(lpmodel, rhs=0, constraints=10)
set.constr.value(lpmodel, rhs=0, constraints=11)
set.constr.type(lpmodel,rep(“=”,7),5:11)
lp.control(lpmodel,sense=’min’)
write.lp(lpmodel,’model.lp’,type=’lp’)
#solve the model, if this return 0 an optimal solution is found
solve(lpmodel)
#this return the proposed solution
get.objective(lpmodel)
get.variables(lpmodel)
実行結果はこんな感じです。
> get.objective(lpmodel)
[1] 301250
> get.variables(lpmodel)
[1] 20000 10000 0 0 0 0 0 15000 15000 0 0 0 0 0 0 0 0 0 0 20000 0 0 10000 0 0
[26] 0 20000 0 0 0 0 10000 0 0 20000 0 0 0
4. 投資問題
生産計画の問題がわかればおそらく単純な投資問題が解けるのはわかるんじゃないかなと思います。
例えば各投資ごとにリスクとリターンが解っている時に、期待収益を最大化させる組み合わせを選ぶとかって感じですね。
なのでここでは時系列を加味した投資問題のモデリングを行います。
とある学生は20000ドルの資金を所持しており、固定された金利の投資の選択を３つ所持しています。
彼のゴールは卒業する3年後に所持している現金の量を最大化する事です。また、彼は中の親しい友人と毎年年末にスキー旅行にいく為にそこで5000ドルを必要とします。
絵にするとこんな感じです。
丸が投資を決めるポイントで、半年ごとに決定する機会がやってきます。（半年ごとなのは投資商品の最短期間が半年だからで、仮に3ヶ月の商品があれば3ヶ月になります）
3年後の現金の保有を最大化したいので、３年後にはすべての投資が帰ってくる状態でなくてはなりません。
さて、これをモデルに書いてみると以下のようになります。
慣れていないと目的関数をどう書くか？でちょっと悩んでしまいますね。
目的が3年後の現金の保有の最大化なので、3年後に現金を生む物を最大化させればよいです。つまり、1期前からの現金の持越し、１期前の6-monthsへの投資量に収益率を掛けたもの、2期前の12-monthsへの投資量に収益率を掛けたもの、4期前の24-monthsへの投資量に収益率を掛けたものの合計からスキー旅行代の5000を差し引いた値になります
あとは1期毎に制約条件を書いてあげるだけです。
1期目の投資と現金持越し額の合計は、保持している資金20000を超えてはいけないというのが1つ目の制約です。
2期目から投資する額の総計は１期目から持ち越された現金と1期目に投資した6-monthsのリターンの合計額以上になってはいけないという制約が必要です。
3期目では先ほどと同じように1期前からの現金と投資に加え、2期前からの12-monthsのリターンも加わります。しかし、ここでスキー旅行に行くので、5000ドルを差し引かねばなりません。よって、入ってきた金額と投資する金額の差が5000になるように、右辺に-5000を置きます。
というのを6期までやってあげると上の様な制約条件が書けるわけです。
で、Rでやってみるとこんな感じになります。
lpmodel<-make.lp(6,20)
#x 1-6
#y 1-5
#w 1-3
#s 1-6
set.column(lpmodel,1,c(1,-1.049), indices=c(1,2))
set.column(lpmodel,2,c(1,-1.049), indices=c(2,3))
set.column(lpmodel,3,c(1,-1.049), indices=c(3,4))
set.column(lpmodel,4,c(1,-1.049), indices=c(4,5))
set.column(lpmodel,5,c(1,-1.049), indices=c(5,6))
set.column(lpmodel,6,c(1,1.049), indices=c(6,0))
set.column(lpmodel,7,c(1,-1.12), indices=c(1,3))
set.column(lpmodel,8,c(1,-1.12), indices=c(2,4))
set.column(lpmodel,9,c(1,-1.12), indices=c(3,5))
set.column(lpmodel,10,c(1,-1.12), indices=c(4,6))
set.column(lpmodel,11,c(1,1.12), indices=c(5,0))
set.column(lpmodel,12,c(1,-1.277), indices=c(1,5))
set.column(lpmodel,13,c(1,-1.277), indices=c(2,6))
set.column(lpmodel,14,c(1,1.277), indices=c(3,0))
set.column(lpmodel,15,c(1,1), indices=c(1,2))
set.column(lpmodel,16,c(1,1), indices=c(2,3))
set.column(lpmodel,17,c(1,1), indices=c(3,4))
set.column(lpmodel,18,c(1,1), indices=c(4,5))
set.column(lpmodel,19,c(1,1), indices=c(5,6))
set.column(lpmodel,20,c(1,1), indices=c(0,1))
set.constr.value(lpmodel, rhs=20000, constraints=1)
set.constr.value(lpmodel, rhs=0, constraints=2)
set.constr.value(lpmodel, rhs=-5000, constraints=3)
set.constr.value(lpmodel, rhs=0, constraints=4)
set.constr.value(lpmodel, rhs=-5000, constraints=5)
set.constr.value(lpmodel, rhs=0, constraints=6)
set.constr.type(lpmodel,
types = rep(3,6),
constraints = 1:6)
lp.control(lpmodel,sense=’max’)
write.lp(lpmodel,’model.lp’,type=’lp’)
#solve the model, if this return 0 an optimal solution is found
solve(lpmodel)
#this return the proposed solution
get.objective(lpmodel)
get.variables(lpmodel)
実行結果はこんな感じです。
実際の解はこの値から５０００を引いたものになります。
> get.objective(lpmodel)
[1] 16619.8
> get.variables(lpmodel)
[1] 0.000 0.000 0.000 0.000 0.000 0.000 4464.286 0.000 0.000 0.000 14839.107 15535.714 0.000 0.000 0.000
[16] 0.000 0.000 0.000 0.000 0.000
12-monthsをひたすら買うって感じですね。
5. Loading of tankers(若干整数計画)
ある石油会社はある国の複数種類の石油商品の需要を満たす長期契約を結んでおり、そのためにタンカーを使って製油所から精油を積み込んでその国に送り込もうとしています。
タンカーには5つのCargo(貯蔵庫)があり、それぞれの許容量が決まっています。
製油所では(91-Octane, 95-Octane, 96-Octane)の3つの商品を生成しており、それらは同じCargoに入れてしまうと混ざってしまう為に別々のCargoに入れる必要があります。
それぞれの商品には需要量がありますが、必ずしもこれを満たす必要はなく、許容値分までなら需要を下回る事が出来ます。しかしこの場合には下回った分のペナルティーコストを相手国に支払う必要が出てきてしまいます。（表のペナルティコストは1単位当たりのコスト）
この会社がペナルティコストをもっとも小さくしようと考えた時、どのような割り振りを行うべきなのでしょうか？
というのがここでのお題です。
この問題ではCargoに複数の商品を混ぜる事が出来ないので、あるCargoに一つの商品を割り振る事が必要になります。
よって整数計画法の考えをここに持ち込みます。
モデルにしてみると下の様になります。
そして、これをRでやるとこんな感じになります。
library(lpSolveAPI)
#used for result visualization
library(ggplot2)
library(reshape)
library(gridExtra)
cargo <- data.frame(
row.names = c(1:5),
capacity = c(5400, 5600, 2200, 1800, 3400)
)
product <- data.frame(
row.names = c(91,95,96),
demand = c(4500,6000,6800),
nondeliver.max =c(800,1800, 1000),
penalty = c(0.5, 0.4, 0.3)
)
lpmodel<-make.lp(26,33)
#y[i,j] 3*5 = 15 : 1 if product i is in cargo j, otherwise 0.
#x[i,j] 3*5 = 15 :volume of product i in cargo j.
#s[i] 3 :volume of product i which is not delivered.
column <- 0
const <- 0
for(j in 1:5){
for(i in 1:3){
column <- column + 1
const <- const + 1
set.column(lpmodel, column, c(1, -1*cargo$capacity[j]), indices = c(j, (const+8)) )
set.type(lpmodel, column, type = “binary”)
}
}
const <- 0
for(j in 1:5){
for(i in 1:3){
column <- column + 1
const <- const + 1
set.column(lpmodel, column, c(1,1), indices = c(i+5, (const + 8)))
}
}
const <- 0
for(i in 1:3){
column <- column + 1
const <- const + 1
set.column(lpmodel, column, c(product$penalty[i], 1, 1), indices = c(0, i+5, i+23))
}
set.constr.value(lpmodel, rhs=rep(1,5), constraints=1:5)
set.constr.value(lpmodel, rhs=product$demand, constraints=6:8)
set.constr.value(lpmodel, rhs=rep(0,15), constraints=9:23)
set.constr.value(lpmodel, rhs=product$nondeliver.max, constraints=24:26)
set.constr.type(lpmodel,
types = rep(3,5),
constraints = 1:5)
set.constr.type(lpmodel,
types = rep(3,3),
constraints = 6:8)
set.constr.type(lpmodel,
types = rep(1,18),
constraints = 9:26)
lp.control(lpmodel,sense=’min’)
write.lp(lpmodel,’model.lp’,type=’lp’)
#solve the model, if this return 0 an optimal solution is found
solve(lpmodel)
#this return the proposed solution
get.objective(lpmodel)
この問題では少しモデルのコーディングの仕方を変えています。毎回数字打ち込んであたふたするのは結構非効率なので。
こちらの方がより大きなモデルを組みたいときに便利です。
まず最初にバイナリ変数のy[i,j]を定義します。この変数は、どの商品がどのカーゴに入るか？を0と1で表してくれます。なのでR上でもset.type()を用いてbinaryである事を定義しています。
column <- 0
const <- 0
for(j in 1:5){
for(i in 1:3){
column <- column + 1
const <- const + 1
set.column(lpmodel, column, c(1, -1*cargo$capacity[j]), indices = c(j, (const+8)) )
set.type(lpmodel, column, type = “binary”)
}
}
次にx[i,j]を定義します。これはCargo[j]の中に入っている商品iの量を表す変数です。
const <- 0
for(j in 1:5){
for(i in 1:3){
column <- column + 1
const <- const + 1
set.column(lpmodel, column, c(1,1), indices = c(i+5, (const + 8)))
}
}
最後に需要を満たさなかった量である、s[i]を定義します。
const <- 0
for(i in 1:3){
column <- column + 1
const <- const + 1
set.column(lpmodel, column, c(product$penalty[i], 1, 1), indices = c(0, i+5, i+23))
}
制約条件も種類ごとにまとめて作っています。
あるカーゴに複数の商品が入らないようにするための制約
set.constr.value(lpmodel, rhs=rep(1,5), constraints=1:5)
配送量と需要を満たせなかった量の合計が、需要量になる為の制約
set.constr.value(lpmodel, rhs=product$demand, constraints=6:8)
ある商品のあるCargoにおける配送量が、需要量×yを上回らないようにする制約。
もしある商品があるCargoを使わないという判断になると、「需要量×y」は0になるので配送量も0になる。使う判断になれば「需要量×y」は需要量となる。
set.constr.value(lpmodel, rhs=rep(0,15), constraints=9:23)
満たさない需要量が許容量を超えないようにする制約。
set.constr.value(lpmodel, rhs=product$nondeliver.max, constraints=24:26)
といった感じで４つの種類の制約条件をそれぞれまとめて作っています。
実行結果はこんな感じです。
> get.objective(lpmodel)
[1] 160
> get.variables(lpmodel)
[1] 1 0 0 0 0 1 0 1 0 0 0 1 0 1 0 4500 0 0
[19] 0 0 5000 0 2200 0 0 0 1800 0 3400 0 0 400 0
最初の15個がy[i,j]で、合計すると５なので、カーゴが５個のみで商品のかぶりなく使われていることが解ります。
6. まとめと参考文献
さて、一旦線形計画法の基礎的な問題はおさらいしたとおもいます。
応用的なものも、結局はこういったアイデアの組合せだったりするので、頑張ればどうにかなります。頑張ってくださいｗ
ただどうしても非線形な問題を扱わないといけなかったり、確率的な要素に配慮しないといけなかったりというケースがあるかと思います。そういった場合には素直に非線形計画とかstochastic optimizationとかやりましょう。
network optimizationとかなんか業務で使えそうだなーと良く思ってるんで使いたいっすね。
※今回の参考資料