2つの目的が共存する組合せ問題の最適解を求めて（Python-MIP）

題材となる問題はこちら
“とあるスーパー、過去3ヵ月以内のレシートを合計N万円分集めると500円券をN枚貰える“

以前、numpyとpandasを用いていましたが、ほぼスクラッチでのグリッドサーチを実装してみました。
こちらがその記事です。

しかし、20枚のレシートで5時間弱の処理時間を要してしまい、実用性に欠けていました。
今回は既存のライブラリを用いて、問題を定式化し、実用性の高いものを作ろうと思います。

#環境
OS：Windows10
CPU：AMD Ryzen 9 5900X 12-Core
メモリ：32GB
バージョン：python 3.9
ツール：jupyter Notebook
※以前のグリッドサーチ実装時と同様です。

#問題文の変換
題材となる問題（以下、レシート問題）は、所謂、0-1ナップサック問題の類似問題と言えます。
一般的な0-1ナップサック問題と比較しやすいように、問題文を変換します。

0-1ナップサック問題

合計重量がW以下かつ、合計価値が最大になる最適な品物の組み合わせは？

レシート問題

合計金額がW以上かつ、合計価値（レシートの古さ）が最大で合計金額が最小になる最適なレシートの組み合わせは？

これでこの問題の要点がわかりやすくなりました。要するに、レシートを一切無駄にすることなく、これからもずっと最大限の利益を上げていきたいのです。
そしてタイトルにもあるように、この問題には目的が2つあることがわかりました。しかし、必ずしも合計価値の最大＝合計金額の最小とは限りません。これを解消するため、今回は2つの方法を試すことにしました。そのことは定式化にて説明します。

#問題の定式化
では、先ほどの問題文を定式化していきます。

0-1ナップサック問題

目的関数（最大）　：　Maximize\sum_{i=1}^{n} v_{i} x_{i} \\
　　制約条件　：　\sum_{i=1}^{n} w_{i} x_{i}\leq W \\
x_{i}\in \{0,1\} (i=1,...,N)

レシート問題（方法①）

目的関数（最大）　：　Maximize\sum_{i=1}^{n} v_{i} x_{i} \\
　制約条件①　：　\sum_{i=1}^{n} w_{i} x_{i}\geq W \\
　　制約条件②　：　\sum_{i=1}^{n} w_{i} x_{i}\leq W+a \\
x_{i}\in \{0,1\} (i=1,...,N)

レシート問題（方法②）

目的関数（最大）　：　Minimize\sum_{i=1}^{n} w_{i} x_{i} \\
　制約条件①　：　\sum_{i=1}^{n} w_{i} x_{i}\geq W \\
　　制約条件②　：　\sum_{i=1}^{n} w_{i} x_{i}\leq W+a \\
　制約条件③　：　\sum_{i=1}^{n} v_{i} x_{i}\geq V \\
x_{i}\in \{0,1\} (i=1,...,N)

どちらか片方を制約に変換することで、2つの方法としました。その過程で、妥協金額aを設けています。
また、合計金額Wに対してVが合計価値で、レシート金額wに対してvがレシート価値となっています。
この2つの方法を実装し、比べてみようと思います。

#価値の付与
実装する前に、まずは各レシートの古さを価値に変換します。
500円券の交換日を0、遡って3ヵ月前を1として、レシートの発行日を正規化します。

import time
import datetime
from dateutil.relativedelta import relativedelta

Useday=datetime.datetime.strptime("2021/11/26","%Y/%m/%d")#交換日
term=Useday-(Useday-relativedelta(months=3))#3ヵ月間
#古い順に並べています。
List_days=["2021/08/28","2021/09/05","2021/09/08",
          "2021/09/15","2021/09/15","2021/09/22",
          "2021/09/29","2021/09/29","2021/09/30",
          "2021/10/05","2021/10/08","2021/10/10",
          "2021/10/20","2021/10/25","2021/10/26",
          "2021/10/30","2021/11/08","2021/11/05",
          "2021/11/08","2021/11/25"
         ]

List_old1=[]
List_old2=[]
List_old3=[]
List_day=[]
for day in List_days:
    a=Useday-datetime.datetime.strptime(day,"%Y/%m/%d")
    List_day.append(a.days)
    List_old1.append(a.days/term.days)#価値パターン1
    List_old2.append((a.days/term.days)**2)#価値パターン2
    List_old3.append((a.days/term.days)**5)#価値パターン3

plt.plot(List_day,List_old1,label="List_old1",marker="o", color="red", linestyle="--")
plt.plot(List_day,List_old2,label="List_old2",marker="o", color="blue", linestyle="--")
plt.plot(List_day,List_old3,label="List_old3",marker="o", color="green", linestyle="--")
plt.legend(loc="upper left")
plt.show()

ここで価値の重みを複数パターン作ってみました。どういう違いになるか実装後に見てみることにします。

#実装
それでは、Python-MIPを用いて実装していきます。
そして今回は、[W,W+a]を反復処理することで最適解がどこにあるか見てみることにします。

from mip import Model, xsum, minimize, maximize, BINARY
from itertools import count

#方法1
def Cal1(w,a,v,W):
    start_time=time.process_time()

    r=range(len(w))
    V=0
    for num in count():
        m=Model("knapsack")
        m.infeas_tol=10**-11 #Python-MIPの演算有効桁数を下げます。
        x=[m.add_var(var_type=BINARY) for i in r]
        m.objective=minimize(xsum(w[i]*x[i] for i in r))
        m+=xsum(w[i]*x[i] for i in r)>=W
        m+=xsum(w[i]*x[i] for i in r)<=W+a
        m+=xsum(v[i]*x[i] for i in r)>=V 
        m.optimize()
        
        if None!=x[0].x:
            selected=[i for i in r if x[i].x>=0.99]
            result1=sum([w[i] for i in selected])
            result2=[w[i] for i in selected]
            result3=[w.index(List_money[i]) for i in selected]
            result4=sum([v[i] for i in selected])
            print("合計金額:{}\nレシート金額:{}\nレシート番号:{}\n合計価値:{}\n".format(result1,result2,result3,result4))
        else:    
            break
        V=result4+v[-1]#今回は最後のレシートが最も価値が低いので利用します。
        
    end_time=time.process_time()
    p_time=end_time-start_time
    print("処理時間:{}\n解答数:{}".format(p_time,num))

#方法2
def Cal2(w,a,v,W):
    start_time=time.process_time()

    r=range(len(w))
    for num in count():
        if a>=0:
            m=Model("knapsack")
            m.infeas_tol=10**-11 #NIPの最低演算桁を下げる
            x=[m.add_var(var_type=BINARY) for i in r]
            m.objective=maximize(xsum(v[i]*x[i] for i in r))
            m+=xsum(w[i]*x[i] for i in r)>=W
            m+=xsum(w[i]*x[i] for i in r)<=W+a
            m.optimize()
        else:
            break
        if None!=x[0].x:
            selected=[i for i in r if x[i].x>=0.99]
            result1=sum([w[i] for i in selected])
            result2=[w[i] for i in selected]
            result3=[w.index(w[i]) for i in selected]
            result4=sum([v[i] for i in selected])
            print("合計金額:{}\nレシート金額:{}\nレシート番号:{}\n合計価値:{}\n".format(result1,result2,result3,result4))
        else:
            break
        a=result1-W-1
        
    end_time=time.process_time()
    p_time=end_time-start_time
    print("処理時間:{}\n解答数:{}".format(p_time,num))

#実行（結果の確認）
それでは実装してみます。レシートは前回と同じなので、同一の解があるはずです。

List_money=[347,2599,1706,1560,2405,5638,
            2876,3375,2261,2746,2579,2067,
            2296,1239,534,1677,1985,989,
            4565,562
           ]
money=40000
a=50
Cal1(List_money,a,List_old1,money) #価値パターン1
Cal2(List_money,a,List_old1,money) 

結果(方法①)

合計金額:40004
レシート金額:[347, 2599, 1560, 2405, 5638, 2876, 3375, 2261, 2746, 2579, 2067, 1239, 534, 1677, 1985, 989, 4565, 562]
レシート番号:[0, 1, 3, 4, 5, 6, 7, 8, 9, 10, 11, 13, 14, 15, 16, 17, 18, 19]
合計価値:9.282608695652172

合計金額:40021
レシート金額:[347, 2599, 1706, 1560, 2405, 5638, 2876, 3375, 2261, 2579, 2067, 2296, 534, 1677, 1985, 989, 4565, 562]
レシート番号:[0, 1, 2, 3, 4, 5, 6, 7, 8, 10, 11, 12, 14, 15, 16, 17, 18, 19]
合計価値:9.630434782608694

合計金額:40033
レシート金額:[347, 2599, 1706, 1560, 2405, 5638, 2876, 3375, 2261, 2746, 2579, 2067, 1239, 534, 1985, 989, 4565, 562]
レシート番号:[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 13, 14, 16, 17, 18, 19]
合計価値:9.84782608695652

処理時間:0.265625
解答数:3

結果(方法②)

合計金額:40033
レシート金額:[347, 2599, 1706, 1560, 2405, 5638, 2876, 3375, 2261, 2746, 2579, 2067, 1239, 534, 1985, 989, 4565, 562]
レシート番号:[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 13, 14, 16, 17, 18, 19]
合計価値:9.84782608695652

合計金額:40021
レシート金額:[347, 2599, 1706, 1560, 2405, 5638, 2876, 3375, 2261, 2579, 2067, 2296, 534, 1677, 1985, 989, 4565, 562]
レシート番号:[0, 1, 2, 3, 4, 5, 6, 7, 8, 10, 11, 12, 14, 15, 16, 17, 18, 19]
合計価値:9.630434782608694

合計金額:40004
レシート金額:[347, 2599, 1560, 2405, 5638, 2876, 3375, 2261, 2746, 2579, 2067, 1239, 534, 1677, 1985, 989, 4565, 562]
レシート番号:[0, 1, 3, 4, 5, 6, 7, 8, 9, 10, 11, 13, 14, 15, 16, 17, 18, 19]
合計価値:9.282608695652172

処理時間:0.296875
解答数:3

1秒未満に処理時間を短縮することができました。先人の英知に学ぶことの重要性を痛感します。
結果としては、どちらの方法でも同じ結果が得られており、グリッドサーチの結果とも同じです。
しかし、目標金額が40000円だとほぼ全てのレシートを使うことになるため、結果の順番が違うだけで、各手法の特性がみれません。
そこで、レシートの金額と目標金額を変更して再度実行してみます。
#実行（手法の比較）
変更点
・1枚目のレシートを347円から5347円に変更
・目標金額を2万円に変更
・妥協金額を10円に変更

List_money=[5347,2599,1706,1560,2405,5638,
            2876,3375,2261,2746,2579,2067,
            2296,1239,534,1677,1985,989,
            4565,562
           ]
money=20000
a=10
Cal1(List_money,a,List_old1,money) #価値パターン1
Cal2(List_money,a,List_old1,money) 

※ここから結果の一部を省略しています。

結果(方法①)

合計金額:20001　レシート番号:[0, 1, 4, 5, 13, 14, 15, 19]
合計金額:20000　レシート番号:[2, 4, 7, 8, 11, 12, 13, 15, 16, 17]
合計金額:20000　レシート番号:[1, 3, 4, 6, 8, 9, 10, 16, 17]
合計金額:20001　レシート番号:[2, 3, 4, 6, 8, 9, 10, 11, 13, 19]
合計金額:20001　レシート番号:[0, 1, 2, 3, 4, 6, 14, 16, 17]
合計金額:20000　レシート番号:[1, 2, 3, 4, 6, 8, 9, 12, 17, 19]
合計金額:20005　レシート番号:[1, 2, 3, 4, 6, 10, 11, 13, 16, 17]
合計金額:20008　レシート番号:[1, 2, 3, 4, 8, 9, 10, 11, 14, 17, 19]
解答数:8

結果(方法②)

合計金額:20008　レシート番号:[1, 2, 3, 4, 8, 9, 10, 11, 14, 17, 19]
合計金額:20005　レシート番号:[1, 2, 3, 4, 6, 10, 11, 13, 16, 17]
合計金額:20000　レシート番号:[1, 2, 3, 4, 6, 8, 9, 12, 17, 19]
解答数:3

方法②で結果が少ないのは、合計金額が同じだと、合計価値が最大の組合せのみを解とするためです。そのため、より最適な組み合わせが知りたければ方法①を使うべきなのでしょうが、合計価値が同じになる確率の方が低いため、処理回数が多くなってしまいます。
そこで、方法②でも最適な組み合わせを得られるようにする工夫が価値の重みの変更です。

#実行（価値パターンの比較）
それでは最後に、方法②を用いて価値パターンについて比較してみます。

List_money=[5347,2599,1706,1560,2405,5638,
            2876,3375,2261,2746,2579,2067,
            2296,1239,534,1677,1985,989,
            4565,562
           ]
money=20000
a=10
Cal2(List_money,a,List_old1,money) #価値パターン1
Cal2(List_money,a,List_old2,money) #価値パターン2
Cal2(List_money,a,List_old3,money) #価値パターン3

結果(方法②)

#価値パターン1
合計金額:20008　レシート番号:[1, 2, 3, 4, 8, 9, 10, 11, 14, 17, 19]
合計金額:20005　レシート番号:[1, 2, 3, 4, 6, 10, 11, 13, 16, 17]
合計金額:20000　レシート番号:[1, 2, 3, 4, 6, 8, 9, 12, 17, 19]
解答数:3

#価値パターン2
合計金額:20008　レシート番号:[0, 1, 2, 3, 4, 8, 10, 17, 19]
合計金額:20001　レシート番号:[0, 1, 2, 3, 4, 6, 14, 16, 17]
合計金額:20000　レシート番号:[1, 2, 3, 4, 6, 8, 9, 12, 17, 19]
解答数:3

#価値パターン3
合計金額:20008　レシート番号:[0, 1, 2, 3, 4, 8, 10, 17, 19]
合計金額:20001　レシート番号:[0, 1, 2, 3, 4, 6, 14, 16, 17]
合計金額:20000　レシート番号:[0, 1, 3, 4, 6, 15, 16, 17, 19]
解答数:3

このように、価値パターンを変えることで、古い方をより優先して選ぶことに成功しました。

#考察とまとめ
今回、目的が2つある組み合わせ問題でしたが、2つの方法に分け、さらに価値の重みも変えて比較してみたところ、方法②と価値パターン3の組合せが最も優れている結果となりました。これでなんとか実用的な運用ができそうです。
最後まで読んで下さり有難うございました。

#参考文献
Python-MIP：https://python-mip.readthedocs.io/en/latest/classes.html
日付や時刻の値への加算・減算と値の比較：https://www.javadrive.jp/python/date/index8.html
Pythonで翌日や翌月みたいな日付の計算をする：https://qiita.com/dkugi/items/8c32cc481b365c277ec2
組合せ最適化を使おう：https://qiita.com/SaitoTsutomu/items/bfbf4c185ed7004b5721

#Webアプリ化
FlaskとRenderを用いてアプリ化しました。

※読み込みには1分程要します。