ナップサック問題入門¶
1. はじめに¶
本資料では、ナップサック問題を例に、CSV ファイルからのデータ読み込み、pulp による定式化と求解、LP ファイルの確認までを一続きで学ぶ。進めるにあたっては以下の点を守ること。
分からないことがあれば、まず検索して調べること。
まずは配られたプログラムや自分が作ったプログラムを全行理解するまで読む（そのために全行コメントを入れている）。
思い通りに動かなければ、関連する変数の中身をprint して確認する。
エラーが出るのであればエラーメッセージを確認する。
エラーメッセージの意味が分からなければコピペして調べる。
知らない関数やメソッド（今回だとstrip やsplit）があったら自分で調べる。
2. CSV ファイルの読み込み¶
問題作成に必要な情報はcsv ファイルに書き込む。 csv ファイルとはカンマで区切られる表形式のファイルである。 表形式なのでExcel で開けるが、諸事情によりExcelでは開かないこと。
この形式のファイルを読み込むところから始める。
# csv ファイル読み込み
id = [] # 品物ID
value = [] # 品物価値
weight = [] # 品物重さ
capacity = 0 # 容器の容量
# with open はファイルを開くおまじない
with open("items.csv", "r") as f: # items.csv をr (read)モードで読み込み
# f にはファイルの中身に関する情報が格納されている
# f.readlines() で行ごとにテキストを取得する
# rows は行を要素とするリストと考えていい
rows = f.readlines() # ファイルの全行をrows に格納
for i in range(len(rows)): # 行数分だけfor 文を回す
if i == 0: continue # 1行目はヘッダー行なのでスキップ
row = rows[i].strip() # strip で各文字列の改行や空白を削除
row = row.split(',') # カンマ区切りで行をリスト化
# "1,2,3".split(",") -> [1, 2, 3]
print(row)
id.append(int(row[0])) # 1 列目はID、整数に変換する（int）
value.append(float(row[1])) # 2 列目は価値、数値に変換する（float）
weight.append(float(row[2]))# 3 列目は重さ、数値に変換する（float）
with open("box.csv", "r") as f: # box.csv をr (read)モードで読み込み
rows = f.readlines() # ファイルの全行をrows に格納
for i in range(len(rows)): # 行数分だけfor 文を回す
if i == 0: continue # 1行目はヘッダー行なのでスキップ
row = rows[i].strip() # strip で各文字列の改行や空白を削除
row = row.split(',') # カンマ区切りで行をリスト化
capacity = int(row[0])
print(id, value, weight, capacity)
['1', '120', '10'] ['2', '250', '21'] ['3', '185', '16'] ['4', '100', '9'] ['5', '130', '12'] ['6', '80', '7'] [1, 2, 3, 4, 5, 6] [120.0, 250.0, 185.0, 100.0, 130.0, 80.0] [10.0, 21.0, 16.0, 9.0, 12.0, 7.0] 60
3. LP ファイルの作成¶
csv ファイルから問題作成に必要なデータを読み込んだら、実際に問題を作成する。 数理最適化問題はLP ファイルと呼ばれる形式で書くことで、ソルバーに読み込ませることができる。 LP ファイルは「ファイル名.lp」という形式で保存する。 .lp というように拡張子はなっているが、中身はただのテキストデータである。 そのため容易にプログラムで生成することができる。 しかし、いちいちテキストファイルを作成するのは大変なので、便利なライブラリを使う。
3.1 pulp による数理最適化問題の定式化¶
pulp とは数理最適化問題を定式化したり、定式化した問題を解くための便利なライブラリである。 以下ではpulp を使って数理最適化問題を定式化する。
- 注意事項
以下、プログラムで使う変数は変数、数理モデルに登場する変数は決定変数と呼ぶ。 混同しないように注意すること。
プログラムと数理モデルでは数の始まりが違う。 プログラムは0 始まり、数理モデルは1 始まりである。 そのためプログラム中の配列は0 始まりで考えるが、lp ファイルでは1 始まりで考える。
まずはpulp を使えるようにするため、pulp をインポートする。
import pulp
ここから定式化を進めていくが、Pulp で定式化する際には、数式で定義した数理モデルとプログラム内の定式化部分で同じ記号を使うようにする、ということに注意すること。
ナップサック問題では定数や記号を以下のように準備した。
- $i$: 品物ID
- $n$: 品物の個数
- $a_i$: 品物𝑖 の重さ [kg]
- $c_i$: 品物𝑖 の輸送価値
- $b$ : 箱の容量 [kg]
これらの記号に合わせてプログラム中の変数を定義しなおす。
n = len(id) # 品物の数（id の長さ）
a_i = weight # 品物i の重さ
c_i = value # 品物i の輸送価値
b = capacity # 箱の容量
名前を変えて定義しなおしただけだが、数理モデルとの対応を取ることで管理が容易になる。
ではここから定式化していく。
3.2 問題の定義¶
ここでは定式化した問題を格納する（プログラムにおける）変数knapsack を定義する。
knapsack に問題の名前と、目的関数を最大化するか最小化するかを設定する。
（目的関数を最大化するか、最小化するか、という設定だけを行っている）
# 問題名（knapsack）と目的関数（最大化か最小化か）の設定
knapsack = pulp.LpProblem('knapsack', pulp.LpMaximize)
pulp.LpProblem とはpulp が持っているLpProblem という名前のメソッドを読んでいる。
メソッドは関数と考えても良い。
なにか入力があると特定の処理をして出力を返してくれる。
例えば、$f(x)=x^2$ という関数は$x=3$ という入力を与えると、$9$ という出力を返してくれる関数である。
python でこの関数を実装すると
def f(x):
return x*x
と書ける。
厳密には違うが、pulp.LpProblem も同じようなものである。
pulp.LpProblem の場合、入力は() 内の'knapsack' とpulp.LpMaximize である（これを引数という）。
これらの入力を与えると、この入力を使って問題を定義して、問題のオブジェクト（問題に関する様々な情報が含まれるもの）が出力として返ってくる。
'knapsack' は問題の名前、pulp.LpMaximize は問題が最大化であることを伝える入力である。
つまり、変数knapsack には'knapsack' という名前の最大化問題オブジェクトが格納される。
以後、このknapsack に目的関数の式、制約条件の式を追加していく。
3.3 変数の定義¶
上で定義した問題knapsack で使う決定変数を定義していく。
決定変数はknapsack とは独立した変数で定義する。
x = [] # 決定変数を格納する変数
for i in range(n):
# LpVariable("決定変数の名前", cat="変数の種類")
x.append(pulp.LpVariable(f'x_{i+1}', cat=pulp.LpBinary))
x
[x_1, x_2, x_3, x_4, x_5, x_6]
例えばx_1 ~x_10 の決定変数があったとしたら、それをもとに決定変数の名前にするとよい。
ただしプログラムは0 始まりなので i+1 としている。
ここで使っているLpBinary は0, 1（バイナリ）しかとらない整数決定変数のことである。
LpBinary 以外にも以下のような変数がある。
LpBinary: バイナリ決定変数LpInteger: 整数値を取る決定変数LpContinuous: 連続値を取る決定
問題に応じて使い分けること。
3.4 目的関数の定義¶
$$ \sum_{i=1}^{n} c_i x_i \rightarrow \max. $$
この目的関数を定義してknapsack に格納する。
# += は 左辺の変数knapsack に右辺の目的関数を定義しろという命令
knapsack += pulp.lpSum(c_i[i] * x[i] for i in range(n))
+= は左辺の問題を格納するための変数に、右辺で定義される目的関数を格納するための命令記号である。
lpSum は$\Sigma$ にあたる総和を取るための、pulp 専用の命令である。
右辺を見ると、変数を数理モデルに合わせて定義しなおしたことで、対応が分かりやすくなっている。
なお、range(n) は0~n-1である。
数理モデルでは$i=1, \dots, n$ だが、プログラムでは配列として$x$ や$c_i$ を扱っているので0~n-1 になる。
ちなみにlpSum の中身は内包表記と呼ばれる書き方になっている。
[ 内包表記 ]
主にfor 分でリストを定義したいときに使ったりする。
0~9 までの値を格納したリストはfor文で
a = []
for i in range(10):
a.append(i)
とすると定義できるが、内包表記をつかうと
a = [i for i in range(10)]
と一文で定義できる。
これを踏まえるとlpSum のかっこの中は内包表記になっていることが分かると思う。
3.5 制約条件の定義¶
$$ \sum_{i=1}^{n} a_i x_i \leq b $$
この制約条件を定義してknapsack に格納する。
knapsack += pulp.lpSum(a_i[i] * x[i] for i in range(n)) <= b
これはknapsack に「重さの合計 ≤ 容量」という制約条件を追加（+=）しろという処理になっている。
ここでも記号を数理モデルに合わせたおかげで対応関係が分かりやすくなっている。
こちらのlpSum の中身も内包表記になっている。
目的関数の定義と違うのは、右辺が式そのもの（xxx <= b のような形）になっている点である。
そのうちこの書き方にも慣れると思う。
これで決定変数、目的関数、制約条件の定義が完了した。 あとは解くだけである。
3.6 ソルバーを利用した数理最適化問題の求解¶
ここまでに定義したknapsack は以下の一文で最適化を実行できる。
早速実行してみる。
# 問題変数.solve() でその問題を解いてくれる
knapsack.solve(pulp.PULP_CBC_CMD(msg=True))
Welcome to the CBC MILP Solver Version: 2.10.10 Build Date: Sep 26 2023 command line - cbc /tmp/4af276de13f64df6acb71db04d1928df-pulp.mps -max -timeMode elapsed -branch -printingOptions all -solution /tmp/4af276de13f64df6acb71db04d1928df-pulp.sol (default strategy 1) At line 2 NAME MODEL At line 3 ROWS At line 6 COLUMNS At line 31 RHS At line 33 BOUNDS At line 40 ENDATA Problem MODEL has 1 rows, 6 columns and 6 elements Coin0008I MODEL read with 0 errors Option for timeMode changed from cpu to elapsed Continuous objective value is 701.667 - 0.03 seconds Cgl0004I processed model has 1 rows, 6 columns (6 integer (6 of which binary)) and 6 elements Cutoff increment increased from 1e-05 to 4.9999 Cbc0038I Initial state - 1 integers unsatisfied sum - 0.333333 Cbc0038I Pass 1: suminf. 0.14286 (1) obj. 699.286 iterations 1 Cbc0038I Solution found of 485 Cbc0038I Rounding solution of 615 is better than previous of 485 Cbc0038I Before mini branch and bound, 4 integers at bound fixed and 0 continuous Cbc0038I Full problem 1 rows 6 columns, reduced to 1 rows 2 columns Cbc0038I Mini branch and bound improved solution from 615 to 635 (0.12 seconds) Cbc0038I Round again with cutoff of 646.167 Cbc0038I Pass 2: suminf. 0.14286 (1) obj. 699.286 iterations 0 Cbc0038I Pass 3: suminf. 0.35533 (1) obj. 646.167 iterations 1 Cbc0038I Pass 4: suminf. 0.03533 (1) obj. 646.167 iterations 1 Cbc0038I Solution found of 655 Cbc0038I Before mini branch and bound, 3 integers at bound fixed and 0 continuous Cbc0038I Full problem 1 rows 6 columns, reduced to 1 rows 3 columns Cbc0038I Mini branch and bound did not improve solution (0.13 seconds) Cbc0038I Round again with cutoff of 668.333 Cbc0038I Pass 5: suminf. 0.14286 (1) obj. 699.286 iterations 0 Cbc0038I Pass 6: suminf. 0.26667 (1) obj. 668.333 iterations 1 Cbc0038I Pass 7: suminf. 0.14286 (1) obj. 699.286 iterations 1 Cbc0038I Pass 8: suminf. 0.26667 (1) obj. 668.333 iterations 1 Cbc0038I Pass 9: suminf. 0.14286 (1) obj. 699.286 iterations 1 Cbc0038I Pass 10: suminf. 0.10256 (1) obj. 668.333 iterations 3 Cbc0038I Pass 11: suminf. 0.38095 (1) obj. 689.762 iterations 2 Cbc0038I Pass 12: suminf. 0.46667 (1) obj. 668.333 iterations 1 Cbc0038I Pass 13: suminf. 0.06667 (1) obj. 668.333 iterations 1 Cbc0038I Solution found of 685 Cbc0038I Before mini branch and bound, 2 integers at bound fixed and 0 continuous Cbc0038I Full problem 1 rows 6 columns, reduced to 1 rows 4 columns Cbc0038I Mini branch and bound did not improve solution (0.13 seconds) Cbc0038I Round again with cutoff of 693.5 Cbc0038I Reduced cost fixing fixed 2 variables on major pass 4 Cbc0038I Pass 14: suminf. 0.18750 (1) obj. 700.312 iterations 1 Cbc0038I Pass 15: suminf. 0.22432 (1) obj. 693.5 iterations 1 Cbc0038I Pass 16: suminf. 0.18750 (1) obj. 700.312 iterations 1 Cbc0038I Pass 17: suminf. 0.22432 (1) obj. 693.5 iterations 1 Cbc0038I Pass 18: suminf. 0.45000 (1) obj. 693.5 iterations 2 Cbc0038I Pass 19: suminf. 0.45000 (1) obj. 693.5 iterations 0 Cbc0038I Pass 20: suminf. 0.37500 (1) obj. 695.625 iterations 2 Cbc0038I Pass 21: suminf. 0.38649 (1) obj. 693.5 iterations 1 Cbc0038I Pass 22: suminf. 0.08500 (1) obj. 693.5 iterations 2 Cbc0038I Pass 23: suminf. 0.08500 (1) obj. 693.5 iterations 0 Cbc0038I Pass 24: suminf. 0.52857 (2) obj. 693.5 iterations 3 Cbc0038I Pass 25: suminf. 0.29615 (1) obj. 693.5 iterations 1 Cbc0038I Pass 26: suminf. 0.29615 (1) obj. 693.5 iterations 0 Cbc0038I Pass 27: suminf. 0.29615 (1) obj. 693.5 iterations 0 Cbc0038I Pass 28: suminf. 0.29615 (1) obj. 693.5 iterations 0 Cbc0038I Pass 29: suminf. 0.08500 (1) obj. 693.5 iterations 1 Cbc0038I Pass 30: suminf. 0.08500 (1) obj. 693.5 iterations 0 Cbc0038I Pass 31: suminf. 0.08500 (1) obj. 693.5 iterations 0 Cbc0038I Pass 32: suminf. 0.38649 (1) obj. 693.5 iterations 2 Cbc0038I Pass 33: suminf. 0.37500 (1) obj. 695.625 iterations 1 Cbc0038I Pass 34: suminf. 0.38649 (1) obj. 693.5 iterations 1 Cbc0038I Pass 35: suminf. 0.38649 (1) obj. 693.5 iterations 0 Cbc0038I Pass 36: suminf. 0.52857 (2) obj. 693.5 iterations 3 Cbc0038I Pass 37: suminf. 0.52857 (2) obj. 693.5 iterations 0 Cbc0038I Pass 38: suminf. 0.29615 (1) obj. 693.5 iterations 1 Cbc0038I Pass 39: suminf. 0.29615 (1) obj. 693.5 iterations 0 Cbc0038I Pass 40: suminf. 0.29615 (1) obj. 693.5 iterations 0 Cbc0038I Pass 41: suminf. 0.18750 (1) obj. 700.312 iterations 3 Cbc0038I Pass 42: suminf. 0.22432 (1) obj. 693.5 iterations 1 Cbc0038I Pass 43: suminf. 0.18750 (1) obj. 700.312 iterations 1 Cbc0038I No solution found this major pass Cbc0038I Before mini branch and bound, 2 integers at bound fixed and 0 continuous Cbc0038I Full problem 1 rows 6 columns, reduced to 1 rows 4 columns Cbc0038I Mini branch and bound did not improve solution (0.13 seconds) Cbc0038I After 0.13 seconds - Feasibility pump exiting with objective of 685 - took 0.01 seconds Cbc0012I Integer solution of 685 found by feasibility pump after 0 iterations and 0 nodes (0.13 seconds) Cbc0038I Full problem 1 rows 6 columns, reduced to 1 rows 3 columns Cbc0006I The LP relaxation is infeasible or too expensive Cbc0013I At root node, 0 cuts changed objective from 701.66667 to 701.66667 in 1 passes Cbc0014I Cut generator 0 (Probing) - 1 row cuts average 0.0 elements, 2 column cuts (2 active) in 0.000 seconds - new frequency is 1 Cbc0014I Cut generator 1 (Gomory) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is -100 Cbc0014I Cut generator 2 (Knapsack) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is -100 Cbc0014I Cut generator 3 (Clique) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is -100 Cbc0014I Cut generator 4 (MixedIntegerRounding2) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is -100 Cbc0014I Cut generator 5 (FlowCover) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is -100 Cbc0014I Cut generator 6 (TwoMirCuts) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is -100 Cbc0014I Cut generator 7 (ZeroHalf) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is -100 Cbc0001I Search completed - best objective 685, took 0 iterations and 0 nodes (0.13 seconds) Cbc0035I Maximum depth 0, 1 variables fixed on reduced cost Cuts at root node changed objective from 701.667 to 701.667 Probing was tried 1 times and created 3 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) ZeroHalf was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Result - Optimal solution found Objective value: 685.00000000 Enumerated nodes: 0 Total iterations: 0 Time (CPU seconds): 0.13 Time (Wallclock seconds): 0.13 Option for printingOptions changed from normal to all Total time (CPU seconds): 0.16 (Wallclock seconds): 0.17
1
何やらメッセージが出てきたが、これは問題を解く過程を表示している。
Result - Optimal solution found というメッセージ以降が大切である。
Result - Optimal solution found
Objective value: 685.00000000
Enumerated nodes: 0
Total iterations: 0
Time (CPU seconds): 0.00
Time (Wallclock seconds): 0.00
Option for printingOptions changed from normal to all
Total time (CPU seconds): 0.00 (Wallclock seconds): 0.01
まずOptimal solution found は最適解が見つかったことを示している。
無事に最適化できたことが分かる。
Objective value: 685.00000000
これは最適化された結果、目的関数の値が685 になった、という意味である。
つまり箱には合計価値685 の品物を入れることができ、それが最大だということである。
この行以降は計算時間に関連する項目である。
問題が小さくほとんど時間がかからなかったことが分かる（表示される時間は実行環境により多少変わる）。
3.7 結果の取得¶
無事最適化できたので、実際にどの決定変数が1 になってどの決定変数が0 になったのか確認する。
print('選んだアイテム:')
for i in range(n): # すべての品目をループ
if x[i].varValue == 1: # varValue で 1 かどうか判定
print(' -', id[i], ', 価値:', value[i], ', 重さ:', weight[i]) # 選ばれた品目名を表示
print('合計重量 = ', sum([a_i[i] * x[i].varValue for i in range(n)])) # 箱に入れた品物の重さの合計
print('合計価値 =', knapsack.objective.value()) # 最終的な価値を表示
選んだアイテム: - 1 , 価値: 120.0 , 重さ: 10.0 - 2 , 価値: 250.0 , 重さ: 21.0 - 3 , 価値: 185.0 , 重さ: 16.0 - 5 , 価値: 130.0 , 重さ: 12.0 合計重量 = 59.0 合計価値 = 685.0
実際に制約条件を満たしていることが確認できる。
x[i].varValue は問題を解き終わった後の品物$i$の決定変数$x_i$ の値を取得している。
最適化された後であれば、決定変数$x_i$ の値は定まっている。
これを最適解と呼ぶ。
注意 最適解とは最終的に得られた決定変数の値（今回でいうと$x_i$ の値）のことである。 つまり今回は $x_1 = 1, x_2 = 1, x_3 = 1, x_4 = 0, x_5 = 1, x_6 = 0$ が最適解である。 目的関数の値（今回だと$685$） は最適値と呼ぶ。 最適解が分かれば、その最適解を目的関数に代入すれば最適値が分かる。
3.8 LP ファイルの取得¶
最後に定式化した問題をLP ファイルに書き出してみる。 以下の一文で出力できる。
knapsack.writeLP('knapsack.lp') # ファイル名.lp
[x_1, x_2, x_3, x_4, x_5, x_6]
実際にLP ファイルの中身を見てみる。
\* knapsack *\
Maximize
OBJ: 120 x_1 + 250 x_2 + 185 x_3 + 100 x_4 + 130 x_5 + 80 x_6
Subject To
_C1: 10 x_1 + 21 x_2 + 16 x_3 + 9 x_4 + 12 x_5 + 7 x_6 <= 60
Binaries
x_1
x_2
x_3
x_4
x_5
x_6
End
OBJ や_C1 はそのあとに続く目的関数や制約式の名前である。
あまり気にしなくていい。
見てみると実際にLP ファイルとして出力できている。
ちなみに今回はLP ファイルを経由せずに最適化したが、LP ファイルに出力しておくと、Python 以外の言語でもLP ファイルを読み込むことで最適化ができる。
つまり、pulp はPython 専用の機能だが、LP ファイルに出力することで共通言語として数理モデルを扱えるわけである。
4. 最後に¶
最後にここまでのプログラムをまとめて全体を俯瞰する。
import pulp
n = len(id) # 品物の数（id の長さ）
a_i = weight # 品物i の重さ
c_i = value # 品物i の輸送価値
b = capacity # 箱の容量
# 問題名（knapsack）と目的関数（最大化か最小化か）の設定
prob = pulp.LpProblem('knapsack', pulp.LpMaximize)
x = [] # 決定変数を格納する変数
for i in range(n):
# LpVariable("決定変数の名前", cat="変数の種類")
x.append(pulp.LpVariable(f'x_{i+1}', cat=pulp.LpBinary))
# 目的関数
prob += pulp.lpSum(c_i[i] * x[i] for i in range(n))
# 制約条件
prob += pulp.lpSum(a_i[i] * x[i] for i in range(n)) <= b
# 最適化
prob.solve(pulp.PULP_CBC_CMD(msg=False)) # msg=False にすると長い出力がなくなる
# 得られた解の出力
print('選んだアイテム:')
for i in range(n): # すべての品目をループ
if x[i].varValue == 1: # varValue で 1 かどうか判定
print(' -', id[i], ', 価値:', value[i], ', 重さ:', weight[i]) # 選ばれた品目名を表示
print('合計重量 = ', sum([a_i[i] * x[i].varValue for i in range(n)])) # 箱に入れた品物の重さの合計
print('合計価値 =', prob.objective.value()) # 最終的な価値を表示
# lp ファイル書き出し
prob.writeLP('knapsack.lp') # ファイル名.lp
選んだアイテム: - 1 , 価値: 120.0 , 重さ: 10.0 - 2 , 価値: 250.0 , 重さ: 21.0 - 3 , 価値: 185.0 , 重さ: 16.0 - 5 , 価値: 130.0 , 重さ: 12.0 合計重量 = 59.0 合計価値 = 685.0
[x_1, x_2, x_3, x_4, x_5, x_6]
5. ビンパッキング問題に適用するにあたって¶
上の例題では要素が1次元の変数(x_1, x_2, ...) を扱った。 ビンパッキング問題では1次元の変数(y_1, y_2, ...) に加えて2次元の変数(x_1,1, x_1,2, ...) も扱う。
そこで以下では2次元変数のプログラム上での定義の仕方をサンプルとして示す。 このまま使うこともできるが、意味をよく理解して使うこと。
# 添え字が二つの変数（x_ij）の定義
# 例示として、i=1,2,...,n / j=1,2,...,m / n=3 / m=5 というケースを考える
'''
添え字が二つの変数は行列と考えるとわかりやすい。
例えばx_i,j でi=1~3, j=1~5 なら
1列目から5列目までを並べると次のようになる（1行目がi=1、2行目がi=2、3行目がi=3）。
x_1,1, x_1,2, x_1,3, x_1,4, x_1,5
x_2,1, x_2,2, x_2,3, x_2,4, x_2,5
x_3,1, x_3,2, x_3,3, x_3,4, x_3,5
こう考えるとx_ij のi が行番号、j が列番号と考えることができる。
リストとして考えると
[[x_1,1, x_1,2, x_1,3, x_1,4, x_1,5],
[x_2,1, x_2,2, x_2,3, x_2,4, x_2,5],
[x_3,1, x_3,2, x_3,3, x_3,4, x_3,5]]
という感じで、リストの中にリストが入っている、2次元リストとなる。
これは、LpVariable でも同じような感じで定義する。
以下を参考にすること。
'''
n = 3 # 実際にはcsv ファイルの行数などがこれにあたる
m = 5 # 実際にはcsv ファイルの行数などがこれにあたる
x = [] # x は2次元配列
# x = [["x_11", "x_12", ...], ["x_21", "x_22", ...], ...]]
x = [] # 決定変数を格納する変数
for i in range(n):
row = [] # i 行目のx を一時的に格納しておくための変数
for j in range(m):
# i 行j 列のx, つまりx_i,j を定義する
# いったんrow に定義したx_i,j を格納する（これでi 行目のリストrow が完成）
row.append(pulp.LpVariable(f'x_{i+1}_{j+1}', cat=pulp.LpBinary))
# for j... で作ったi 行目のrow をx に追加する
x.append(row)
x
[[x_1_1, x_1_2, x_1_3, x_1_4, x_1_5], [x_2_1, x_2_2, x_2_3, x_2_4, x_2_5], [x_3_1, x_3_2, x_3_3, x_3_4, x_3_5]]