docplexのサンプルのwarehouseを読み解きます。
以下に英語のマニュアルで解説がありますが、この記事では日本語で解説します。

実行環境
- CPLEX 22.1.1
- docplex-2.29.241
- Windows 11 64bit

問題

既存の店舗に供給するための倉庫を開設する場所をいくつか検討しています。それぞれの倉庫には、対応できる店舗数を示す許容店舗数と開設固定費があります。各店舗は1つの倉庫からしか供給を受けず、店舗への供給コストは選択した倉庫によって異なります。最もコストを下げられるように、どの倉庫を開設し、どの店舗にどの倉庫から供給するかを決定します。

利用データ

各倉庫から供給できる店舗数上限:

	ボン	ボルドー	ブリュッセル	ロンドン	パリ	ローマ
許容店舗数	1	4	2	2	1	2

2.各倉庫の開設固定費

	ボン	ボルドー	ブリュッセル	ロンドン	パリ	ローマ
開設固定費	20	35	28	32	37	25

3.各倉庫から各店舗への輸送費

	ボン	ボルドー	ブリュッセル	ロンドン	パリ	ローマ
店舗1輸送費	20	24	30	11	25	30
店舗2輸送費	28	27	74	82	83	74
店舗3輸送費	74	97	70	71	96	70
店舗4輸送費	2	55	61	73	69	61
店舗5輸送費	46	96	34	59	83	54
店舗6輸送費	42	22	59	29	67	59
店舗7輸送費	1	5	56	73	59	56
店舗8輸送費	10	73	96	13	43	96
店舗9輸送費	93	35	46	63	85	46
店舗10輸送費	47	65	95	55	71	95

目的関数

「開設固定費＋各倉庫から各店舗への輸送費＝総コスト」を最小化

決定変数

各倉庫を開設するかどうか
どの倉庫からどの店舗に供給するか

制約

各店舗は一つの倉庫から必ず供給を受ける
開設した倉庫からのみ供給する
各倉庫から供給できる店舗数には上限がある

問題の種類

整数計画

解

先に解を紹介します。

目的関数

合計費用：433

決定変数

各倉庫を開設するかどうか：青いセルで表示した倉庫
どの倉庫からどの店舗に供給するか:緑のセルの倉庫と店舗の組み合わせ

この解で面白いのは以下のような点になります。

「ボン」は「店舗7」への費用が1で最低だが、「店舗4」を優先した方が合計では低くなる
「店舗1」は「ブリュッセル」よりも「パリ」から供給した方がコストは安いが、「パリ」の倉庫の開設費用を抑えた方がトータル費用は抑えることができる

pythonコード解説

pythonコードの解説を行っていきます。

ライブラリーのインポート

まず、ライブラリーのインポートをして環境を確認します。
docplex.mpでdocplexの数理計画モデリングのライブラリーをインポートします。

ライブラリーのインポート

try:
    import docplex.mp
except:
    !pip install docplex
from docplex.mp.environment import Environment
env = Environment()
env.print_information()

結果例

* system is: Windows 64bit
* Python version 3.10.9, located at: c:\Users\xxxxx\AppData\Local\Programs\Python\Python310\python.exe
* docplex is present, version is 2.29.241
* CPLEX library is present, version is 22.1.1.0, located at: c:\Users\xxxxx\AppData\Local\Programs\Python\Python310\lib\site-packages
* pandas is present, version is 1.5.2

変数定義

倉庫に関するデータ定義を行っています。

WAREHOUSESというディクショナリを定義しています。
key部は倉庫名、許容店舗数、開設固定費のタプルになっていて、value部は各店舗との輸送費のリストになっています。

データ定義

#倉庫情報：許容店舗数、開設固定費、各店舗と倉庫間の輸送費
WAREHOUSES = {
    ('Bonn', 1, 20): [20, 28, 74, 2, 46, 42, 1, 10, 93, 47],
    ('Bordeaux', 4, 35): [24, 27, 97, 55, 96, 22, 5, 73, 35, 65],
    ('Brussels', 2, 28): [30, 74, 70, 61, 34, 59, 56, 96, 46, 95],
    ('London', 2, 32): [11, 82, 71, 73, 59, 29, 73, 13, 63, 55],
    ('Paris', 1, 37): [25, 83, 96, 69, 83, 67, 59, 43, 85, 71],
    ('Rome', 2, 25): [30, 74, 70, 61, 54, 59, 56, 96, 46, 95]
}

print(WAREHOUSES)
print("#candidate warehouses=%d" % len(WAREHOUSES))
print("#stores to supply=%d" % len(list(WAREHOUSES.values())[0]))

結果

{('Bonn', 1, 20): [20, 28, 74, 2, 46, 42, 1, 10, 93, 47], ('Bordeaux', 4, 35): [24, 27, 97, 55, 96, 22, 5, 73, 35, 65], ('Brussels', 2, 28): [30, 74, 70, 61, 34, 59, 56, 96, 46, 95], ('London', 2, 32): [11, 82, 71, 73, 59, 29, 73, 13, 63, 55], ('Paris', 1, 37): [25, 83, 96, 69, 83, 67, 59, 43, 85, 71], ('Rome', 2, 25): [30, 74, 70, 61, 54, 59, 56, 96, 46, 95]}
#candidate warehouses=6
#stores to supply=10

WAREHOUSESをData Wranglerでみると以下のようになっています。

結果

例えば、以下のようなことがわかります。

ロンドンの倉庫は最大 2 店舗で供給可能
ロンドン倉庫の開設コストは32です
ロンドンからの店舗1への供給にかかるコストは 11 です。

なお、list(WAREHOUSES.values())は店舗への輸送費を二次元配列で取り出しています。

list(WAREHOUSES.values())

[[20, 28, 74, 2, 46, 42, 1, 10, 93, 47],
 [24, 27, 97, 55, 96, 22, 5, 73, 35, 65],
 [30, 74, 70, 61, 34, 59, 56, 96, 46, 95],
 [11, 82, 71, 73, 59, 29, 73, 13, 63, 55],
 [25, 83, 96, 69, 83, 67, 59, 43, 85, 71],
 [30, 74, 70, 61, 54, 59, 56, 96, 46, 95]]

ここから[0]であるBonの配列を取り出して、lenでサイズを出して、店舗数を算出しています。

利用データ

次にモデルで利用する形に変数を作っていきます。
名前付きタプルで倉庫と店舗のクラスを作っています。これはあとで出てくる決定変数の定義や制約の定義を行いやすくするためです。
倉庫はidと許容量と開設固定費からなる名前付きタプルです。
店舗はidのみをもつ名前付きタプルです。
いずれもstr化するとidが戻ります。

名前付きタプルでクラス作成

from collections import namedtuple
# 倉庫のクラス
class TWareHouse(namedtuple("TWarehouse1", ["id", "capacity", "fixed_cost"])):
    def __str__(self):
        return self.id
# 店舗IDのクラス
class TStore(namedtuple("TStore1", ["id"])):
    def __str__(self):
        return 'store_%d' % self.id

以下で定数データを定義しています。

店舗数nb_storesは先ほどのlen(list(WAREHOUSES.values())[0])で取得してもよいと思うのですが、Python2互換を意識してlen(next(itervalues(WAREHOUSES))で取得しています。
輸送コストsupply_costsはWAREHOUSESディクショナリをそのまま使用しています。
倉庫warehousesは、上で定義した名前付きタプルにWAREHOUSESディクショナリのKey部のidと許容量と開設固定費を入れています。
店舗Istoresは、上で定義した名前付きタプルに1－10までの整数をいれています。

定数データの定義

from six import itervalues

# 店舗数
nb_stores = 0 if not WAREHOUSES else len(next(itervalues(WAREHOUSES)))
# 倉庫;許容店舗数と開設固定費を含む名前付きタプルとして定義
warehouses = [TWareHouse(*wrow) for wrow in WAREHOUSES.keys()]
# 輸送コスト;倉庫ディクショナリをそのまま使用
supply_costs = WAREHOUSES
# 店舗IDのリスト;便宜上、storesは0からではなく1からNSTORESまでカウントします
stores = [TStore(idx) for idx in range(1, nb_stores + 1)]

warehousesは以下のようなタプルになります。

同様にstoresも以下のようなタプルになります。

TWareHouseもstoresには以下のように__str__が定義されていましたので、

class TWareHouse(namedtuple("TWarehouse1", ["id", "capacity", "fixed_cost"])):
    def __str__(self):
        return self.id
# 店舗IDのクラス
class TStore(namedtuple("TStore1", ["id"])):
    def __str__(self):
        return 'store_%d' % self.id

Printとするとidが返ります。

名前付きタプルの中身確認

print(warehouses[0])
print(stores[0])

結果

Bonn
store_1

モデルの作成

以下でモデルオブジェクトを作ります。

モデルオブジェクト作成

from docplex.mp.model import Model
warehouse_model = Model()

決定変数

次に決定変数の定義を行っています。

決定変数1: 倉庫の開設の有無

まず、倉庫ごとに 1 つのバイナリ (はい/いいえ) 決定変数を作成します。この変数は、倉庫が開いている場合にのみ 1 になります。モデルオブジェクトの binary_var_dict メソッドを使用して、名前付きタプルwarehousesからディクショナリを作成します。最初の引数keysは、辞書のキーがwarehousesになることを指定しています。 2 番目の引数は単純な文字列 'open' です。これは決定変数の名前を生成するために使用されます。決定変数の名前は、 'open' ＋'_'+名前付きタプルwarehousesの文字列表現 (例えばBonnなどです。このためにTWareHouseの__str__メソッドを定義していました)になります。

「倉庫の開設の有無」決定変数の定義

open_vars = warehouse_model.binary_var_dict(keys=warehouses, name='open')

以下のようなディクショナリが出来上がります。「open_Bonn」などが0か1が入る決定変数になります。

決定変数2: 倉庫と店舗の供給の有無

次に、ペア (倉庫、店舗) ごとに 1 つのバイナリ変数を定義します。これは、binary_var_matrix メソッドを使用して行われます。この場合、変数へのキーは (warehouses、stores) のペアになります。3 番目の引数は単純な文字列 supplyです。

これは決定変数の名前を生成するために使用されます。決定変数の名前は、 'supply' ＋''+名前付きタプルwarehousesの文字列表現＋''+名前付きタプルstoresの文字列表現になります。

「倉庫と店舗の供給の有無」決定変数の定義

supply_vars = warehouse_model.binary_var_matrix(keys1=warehouses, keys2=stores, name='supply')

以下のようなディクショナリが出来上がります。TWareHouse(id='Bonn', capacity=1, fixed_cost=20), TStore(id=1)などがwarehouses、storesのペアのキーで、supply_Bonn_store_1などが0か1が入る決定変数になります。

定義した決定変数を確認してみます。

決定変数の定義の確認

warehouse_model.print_information()

number of variables: 66で- binary=66となっています。
倉庫が6個と店舗が 10個あるため、6*10 + 6 = 66 個のバイナリ変数が定義されていることがわかります。

結果

Model: docplex_model1
 - number of variables: 66
   - binary=66, integer=0, continuous=0
 - number of constraints: 0
   - linear=0
 - parameters: defaults
 - objective: none
 - problem type is: MILP

制約

次に制約の定義を行っています。

制約1: 各店舗は一つの倉庫

「各店舗は一つの倉庫」の制約は以下のように定義されます。
すべての店舗で、利用する倉庫をかならず一つ割り当てるという意味です。

for s in storesで店舗数分の制約を作っています。
sum(supply_vars[w, s] for w in warehouses)は、ある店舗に供給するすべての倉庫数の合計をだしています。
店舗は必ず一つの倉庫から納入するので「==1」の制約になります。
ctname=f'ct01_{str(s)}'で制約にラベルを付けています。

制約1: 各店舗は一つの倉庫

for s in stores:
    warehouse_model.add_constraint(
        warehouse_model.sum(supply_vars[w, s] for w in warehouses) == 1,
        ctname=f'ct01_{str(s)}')

warehouse_model.print_information()

店舗数分の10個の制約が作られています。

数式で表すと以下です。

\sum_{}^{w}Supply_{sw}=1　 \{すべての店舗sに対して\}\\

sumの部分が展開されてsupplyの各決定変数の足し算になっています。

制約1の表示

for s in stores:
    print(warehouse_model.get_constraint_by_name(f'ct01_{str(s)}'))

制約1の表示

ct01_store_1: supply_Bonn_store_1+supply_Bordeaux_store_1+supply_Brussels_store_1+supply_London_store_1+supply_Paris_store_1+supply_Rome_store_1 == 1
ct01_store_2: supply_Bonn_store_2+supply_Bordeaux_store_2+supply_Brussels_store_2+supply_London_store_2+supply_Paris_store_2+supply_Rome_store_2 == 1
ct01_store_3: supply_Bonn_store_3+supply_Bordeaux_store_3+supply_Brussels_store_3+supply_London_store_3+supply_Paris_store_3+supply_Rome_store_3 == 1
ct01_store_4: supply_Bonn_store_4+supply_Bordeaux_store_4+supply_Brussels_store_4+supply_London_store_4+supply_Paris_store_4+supply_Rome_store_4 == 1
ct01_store_5: supply_Bonn_store_5+supply_Bordeaux_store_5+supply_Brussels_store_5+supply_London_store_5+supply_Paris_store_5+supply_Rome_store_5 == 1
ct01_store_6: supply_Bonn_store_6+supply_Bordeaux_store_6+supply_Brussels_store_6+supply_London_store_6+supply_Paris_store_6+supply_Rome_store_6 == 1
ct01_store_7: supply_Bonn_store_7+supply_Bordeaux_store_7+supply_Brussels_store_7+supply_London_store_7+supply_Paris_store_7+supply_Rome_store_7 == 1
ct01_store_8: supply_Bonn_store_8+supply_Bordeaux_store_8+supply_Brussels_store_8+supply_London_store_8+supply_Paris_store_8+supply_Rome_store_8 == 1
ct01_store_9: supply_Bonn_store_9+supply_Bordeaux_store_9+supply_Brussels_store_9+supply_London_store_9+supply_Paris_store_9+supply_Rome_store_9 == 1
ct01_store_10: supply_Bonn_store_10+supply_Bordeaux_store_10+supply_Brussels_store_10+supply_London_store_10+supply_Paris_store_10+supply_Rome_store_10 == 1

図にすると以下のようなイメージになります。Store3まで表示しています。

制約2: 開設した倉庫のみを利用

「開設した倉庫のみを利用」の制約は以下のように定義されます。

for w in warehouses:とfor s in storesで倉庫数×店舗数分の制約を作っています。
供給の有無を示すsupply_vars[w, s]には「0」か「1」が入ります。開設するかどうかのopen_vars[w]にも「0」か「1」が入ります。開設していない（open_vars[w]=0）場合にはその倉庫からの供給は受けられないので、その倉庫からすべての店舗への供給のsupply_vars[w, s] はすべて0でなければなりません。倉庫が開設される（open_vars[w]=1）場合も、店舗が供給を受けるとは限らないので「＝」ではなく「<=」になります。

f'ct02_{str(w)}_{str(s)}'で制約にラベルを付けています。

制約2: 開設した倉庫のみを利用.mod

for w in warehouses:
    for s in stores:
        warehouse_model.add_constraint(
            supply_vars[w, s] <= open_vars[w],
            ctname=f'ct02_{str(w)}_{str(s)}')

warehouse_model.print_information()

店舗数×倉庫数分で60個の制約が作られています。

数式で表すと以下です。

Supply_{sw} <= Open_w \{すべての倉庫wとすべての店舗sに対して\}\\

制約2の表示

for w in warehouses:
    for s in stores:
        print(warehouse_model.get_constraint_by_name(f'ct02_{str(w)}_{str(s)}'))

制約2の表示

ct02_Bonn_store_1: supply_Bonn_store_1 <= open_Bonn
ct02_Bonn_store_2: supply_Bonn_store_2 <= open_Bonn
ct02_Bonn_store_3: supply_Bonn_store_3 <= open_Bonn
ct02_Bonn_store_4: supply_Bonn_store_4 <= open_Bonn
ct02_Bonn_store_5: supply_Bonn_store_5 <= open_Bonn
ct02_Bonn_store_6: supply_Bonn_store_6 <= open_Bonn
ct02_Bonn_store_7: supply_Bonn_store_7 <= open_Bonn
ct02_Bonn_store_8: supply_Bonn_store_8 <= open_Bonn
ct02_Bonn_store_9: supply_Bonn_store_9 <= open_Bonn
ct02_Bonn_store_10: supply_Bonn_store_10 <= open_Bonn
ct02_Bordeaux_store_1: supply_Bordeaux_store_1 <= open_Bordeaux
ct02_Bordeaux_store_2: supply_Bordeaux_store_2 <= open_Bordeaux
ct02_Bordeaux_store_3: supply_Bordeaux_store_3 <= open_Bordeaux
ct02_Bordeaux_store_4: supply_Bordeaux_store_4 <= open_Bordeaux
ct02_Bordeaux_store_5: supply_Bordeaux_store_5 <= open_Bordeaux
ct02_Bordeaux_store_6: supply_Bordeaux_store_6 <= open_Bordeaux
ct02_Bordeaux_store_7: supply_Bordeaux_store_7 <= open_Bordeaux
ct02_Bordeaux_store_8: supply_Bordeaux_store_8 <= open_Bordeaux
ct02_Bordeaux_store_9: supply_Bordeaux_store_9 <= open_Bordeaux
ct02_Bordeaux_store_10: supply_Bordeaux_store_10 <= open_Bordeaux

-----中略-----

ct02_Rome_store_1: supply_Rome_store_1 <= open_Rome
ct02_Rome_store_2: supply_Rome_store_2 <= open_Rome
ct02_Rome_store_3: supply_Rome_store_3 <= open_Rome
ct02_Rome_store_4: supply_Rome_store_4 <= open_Rome
ct02_Rome_store_5: supply_Rome_store_5 <= open_Rome
ct02_Rome_store_6: supply_Rome_store_6 <= open_Rome
ct02_Rome_store_7: supply_Rome_store_7 <= open_Rome
ct02_Rome_store_8: supply_Rome_store_8 <= open_Rome
ct02_Rome_store_9: supply_Rome_store_9 <= open_Rome
ct02_Rome_store_10: supply_Rome_store_10 <= open_Rome

図にすると以下のようなイメージになります(BordeauxとParisのみ表示)。

制約3: 倉庫許容量以下の供給

「倉庫許容量以下の供給」の制約は以下のように定義されます。各倉庫の許容量以下の供給であることを示しています。

for w in warehouses:で倉庫数分の制約を作っています。

sum(supply_vars[w, s] for s in stores)はある倉庫に供給する店舗すべての供給量です。これがその倉庫の許容量w.capacity以下である必要があります。

制約3: 倉庫許容量以下の供給

for w in warehouses:
    warehouse_model.add_constraint(
        warehouse_model.sum(supply_vars[w, s] for s in stores) <= w.capacity,
        ctname=f'ct03_{str(w)}')

warehouse_model.print_information()

倉庫数分の6個の制約が作られています。

数式で表すと以下です。

\sum_{}^{s}Supply_{sw}<=Capacity_w　 \{すべての倉庫wに対して\}\\

制約3の表示

for w in warehouses:
    print(warehouse_model.get_constraint_by_name(f'ct03_{str(w)}'))

制約3の表示

ct03_Bonn: supply_Bonn_store_1+supply_Bonn_store_2+supply_Bonn_store_3+supply_Bonn_store_4+supply_Bonn_store_5+supply_Bonn_store_6+supply_Bonn_store_7+supply_Bonn_store_8+supply_Bonn_store_9+supply_Bonn_store_10 <= 1
ct03_Bordeaux: supply_Bordeaux_store_1+supply_Bordeaux_store_2+supply_Bordeaux_store_3+supply_Bordeaux_store_4+supply_Bordeaux_store_5+supply_Bordeaux_store_6+supply_Bordeaux_store_7+supply_Bordeaux_store_8+supply_Bordeaux_store_9+supply_Bordeaux_store_10 <= 4
ct03_Brussels: supply_Brussels_store_1+supply_Brussels_store_2+supply_Brussels_store_3+supply_Brussels_store_4+supply_Brussels_store_5+supply_Brussels_store_6+supply_Brussels_store_7+supply_Brussels_store_8+supply_Brussels_store_9+supply_Brussels_store_10 <= 2
ct03_London: supply_London_store_1+supply_London_store_2+supply_London_store_3+supply_London_store_4+supply_London_store_5+supply_London_store_6+supply_London_store_7+supply_London_store_8+supply_London_store_9+supply_London_store_10 <= 2
ct03_Paris: supply_Paris_store_1+supply_Paris_store_2+supply_Paris_store_3+supply_Paris_store_4+supply_Paris_store_5+supply_Paris_store_6+supply_Paris_store_7+supply_Paris_store_8+supply_Paris_store_9+supply_Paris_store_10 <= 1
ct03_Rome: supply_Rome_store_1+supply_Rome_store_2+supply_Rome_store_3+supply_Rome_store_4+supply_Rome_store_5+supply_Rome_store_6+supply_Rome_store_7+supply_Rome_store_8+supply_Rome_store_9+supply_Rome_store_10 <= 2

図にすると以下のようなイメージになります(Bonn、BordeauxとParisのみ表示)。

目的関数

次に目的関数の定義を行っています。

決定式1:開設コスト

まず、決定式で開設コストを定義します。
倉庫の開設open_vars[w]×開設固定費w.fixed_costの合計を出しています。

開設コスト

total_opening_cost = warehouse_model.sum(open_vars[w] * w.fixed_cost for w in warehouses)
print(total_opening_cost)

sumが展開された形になります。

結果

20open_Bonn+35open_Bordeaux+28open_Brussels+32open_London+37open_Paris+25open_Rome

図で表すと以下になります。

決定式2:輸送コスト

次に、決定式で輸送コストを定義します。
決定変数「倉庫と店舗の供給の有無」supply_vars×輸送コストsupply_costsの合計を出しています。輸送コストsupply_costsは('Bonn', 1, 20): [20, 28, 74, 2, 46, 42, 1, 10, 93, 47]というようなディクショナリです。supply_costsのリストの添え字は0から始まるため「店舗のid-1」をしています。
「倉庫と店舗の供給の有無」と「輸送コスト」を総当たりで掛け算した合計ということになります。つまり供給の総コストになります。

輸送コスト

total_supply_cost = warehouse_model.sum(
    [supply_vars[w,s] * supply_costs[w][s.id - 1] for w in warehouses for s in stores])

print(total_supply_cost)

sumが展開された形になります（とても長いです）。

結果

20supply_Bonn_store_1+28supply_Bonn_store_2+74supply_Bonn_store_3+2supply_Bonn_store_4+46supply_Bonn_store_5+42supply_Bonn_store_6+supply_Bonn_store_7+10supply_Bonn_store_8+93supply_Bonn_store_9+47supply_Bonn_store_10+24supply_Bordeaux_store_1+27supply_Bordeaux_store_2+97supply_Bordeaux_store_3+55supply_Bordeaux_store_4+96supply_Bordeaux_store_5+22supply_Bordeaux_store_6+5supply_Bordeaux_store_7+73supply_Bordeaux_store_8+35supply_Bordeaux_store_9+65supply_Bordeaux_store_10+30supply_Brussels_store_1+74supply_Brussels_store_2+70supply_Brussels_store_3+61supply_Brussels_store_4+34supply_Brussels_store_5+59supply_Brussels_store_6+56supply_Brussels_store_7+96supply_Brussels_store_8+46supply_Brussels_store_9+95supply_Brussels_store_10+11supply_London_store_1+82supply_London_store_2+71supply_London_store_3+73supply_London_store_4+59supply_London_store_5+29supply_London_store_6+73supply_London_store_7+13supply_London_store_8+63supply_London_store_9+55supply_London_store_10+25supply_Paris_store_1+83supply_Paris_store_2+96supply_Paris_store_3+69supply_Paris_store_4+83supply_Paris_store_5+67supply_Paris_store_6+59supply_Paris_store_7+43supply_Paris_store_8+85supply_Paris_store_9+71supply_Paris_store_10+30supply_Rome_store_1+74supply_Rome_store_2+70supply_Rome_store_3+61supply_Rome_store_4+54supply_Rome_store_5+59supply_Rome_store_6+56supply_Rome_store_7+96supply_Rome_store_8+46supply_Rome_store_9+95supply_Rome_store_10

図で表すと以下になります。

目的関数

「決定式1:開設コスト」と「決定式2:輸送コスト」の和を最小化します。

minimizeで最小化問題と定義します。

目的関数

warehouse_model.minimize(total_opening_cost + total_supply_cost)
warehouse_model.print_information()

数式で表すと以下です。

minimize \sum_{}^{w} Open_w*Fixed+\sum_{}^{w}\sum_{}^{s}SupplyCost_{sw}*{Supply}_{sw}

モデルを解く

ここまでで作成したモデルをsolveで求解します。
log_output=Trueのオプションでログを出力できます。その他のオプションはマニュアルを参照してください。

求解

sol=warehouse_model.solve(log_output=True)

0.06秒で解くことができました。

結果

Version identifier: 22.1.1.0 | 2022-11-27 | 9160aff4d
CPXPARAM_Read_DataCheck                          1
Tried aggregator 1 time.
MIP Presolve eliminated 20 rows and 0 columns.
MIP Presolve modified 2 coefficients.
Reduced MIP has 56 rows, 66 columns, and 202 nonzeros.
Reduced MIP has 66 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.00 sec. (0.16 ticks)
Found incumbent of value 597.000000 after 0.02 sec. (0.35 ticks)
Probing time = 0.00 sec. (0.05 ticks)
Tried aggregator 1 time.
Detecting symmetries...
Reduced MIP has 56 rows, 66 columns, and 202 nonzeros.
Reduced MIP has 66 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.00 sec. (0.16 ticks)
Probing time = 0.00 sec. (0.05 ticks)
Clique table members: 52.
MIP emphasis: balance optimality and feasibility.
MIP search method: dynamic search.
Parallel mode: deterministic, using up to 8 threads.
Root relaxation solution time = 0.00 sec. (0.14 ticks)

        Nodes                                         Cuts/
   Node  Left     Objective  IInf  Best Integer    Best Bound    ItCnt     Gap

*     0+    0                          597.0000        0.0000           100.00%
*     0+    0                          516.0000        0.0000           100.00%
*     0+    0                          491.0000        0.0000           100.00%
      0     0      412.7500    13      491.0000      412.7500       28   15.94%
*     0+    0                          456.0000      412.7500             9.48%
*     0+    0                          433.0000      412.7500             4.68%
      0     0      425.0000    16      433.0000       Cuts: 7       37    1.85%
      0     0        cutoff            433.0000                     44    0.00%
Elapsed time = 0.06 sec. (1.91 ticks, tree = 0.01 MB, solutions = 5)

Implied bound cuts applied:  1
Zero-half cuts applied:  2
Lift and project cuts applied:  1
Gomory fractional cuts applied:  1

Root node processing (before b&c):
  Real time             =    0.06 sec. (1.92 ticks)
Parallel b&c, 8 threads:
  Real time             =    0.00 sec. (0.00 ticks)
  Sync time (average)   =    0.00 sec.
  Wait time (average)   =    0.00 sec.
                          ------------
Total (root+branch&cut) =    0.06 sec. (1.92 ticks)

解をprint_solutionで表示します。
print_zeros=Trueにすると0の決定変数も表示されます。その他のオプションはマニュアルを参照してください。

解の表示

warehouse_model.print_solution()

目的関数は433です。１になった決定変数も表示されます。

結果

objective: 433
status: OPTIMAL_SOLUTION(2)
  open_Bonn=1
  open_Bordeaux=1
  open_Brussels=1
  open_London=1
  open_Rome=1
  supply_Bonn_store_4=1
  supply_Bordeaux_store_2=1
  supply_Bordeaux_store_6=1
  supply_Bordeaux_store_7=1
  supply_Bordeaux_store_9=1
  supply_Brussels_store_1=1
  supply_Brussels_store_5=1
  supply_London_store_8=1
  supply_London_store_10=1
  supply_Rome_store_3=1

結果取得

結果は出ましたが、表示をするだけではなく、変数に代入して自由に利用できるようにしていきます。

目的関数の結果

目的関数の結果はobjective_valueで取得できます。

目的関数の結果

obj = warehouse_model.objective_value
print("optimal objective is %g" % obj)

結果

optimal objective is 433

決定式の結果

決定式の結果はsolution_valueで取得できます。

後でつかいたいのでresultsに入れておきます。

決定式

opening_cost_obj = total_opening_cost.solution_value
supply_cost_obj = total_supply_cost.solution_value
print("total opening cost=%g" % opening_cost_obj)
print("total supply cost=%g" % supply_cost_obj)

# store results for later on...
results=[]
results.append((opening_cost_obj, supply_cost_obj))

結果

total opening cost=140
total supply cost=293

決定変数1: 倉庫の開設の有無

決定変数の結果は、決定式と同様にsolution_valueで取得できます。
ただ、open_varsはbinary_var_dictなので、辞書形式のvaluesに入っています。
表形式にした方がわかりやすいので、keyの名前付きタプルwarehousesからidも取り出して、solution_valueとあわせて、DataFrameにしています。

決定変数1: 倉庫の開設の有無

import pandas as pd

open_vars_w = [w.id for w in open_vars.keys()]
open_vars_v = [round(v.solution_value) for v in open_vars.values()]

pd.DataFrame({"warehouse": open_vars_w,
              "open": open_vars_v})

決定変数2: 倉庫と店舗の供給の有無

supply_varsはbinary_var_matrixなので、決定変数は辞書形式のvaluesに入っています。
表形式にした方がわかりやすいので、keyの名前付きタプルwarehousesのidとstoresのidも取り出して、solution_valueとあわせて、DataFrameにしています。
binary_var_matrixのkeysは二次元なので、k[0]にはwarehouses、k[1]にはstoresが入っています。

決定変数2: 倉庫と店舗の供給の有無

supply_vars_w = [k[0].id for k in supply_vars.keys()]
supply_vars_s = [k[1].id for k in supply_vars.keys()]
supply_vars_v = [round(v.solution_value) for v in supply_vars.values()]

df_suuply = pd.DataFrame({"warehouse": supply_vars_w,
                          "store": supply_vars_s,
                         "supply": supply_vars_v})
df_suuply

以下のようになります。Bordeauxまでを表示しています。

倉庫と店舗でPivotするとよりわかりやすく見ることができます。

決定変数2: 倉庫と店舗の供給の有無 Pivot

df_suuply.pivot_table(values="supply", index="store", 
                     columns="warehouse")

制約2のスラック表示

「スラック」とは制約式の余裕量を示しています。

制約の式はto_stringで抜き出しています。
_extract_after_colonはコロン以降の文字列を抜き出すローカル関数です。

get_constraint_by_name(f'ct02_{str(w)}_{str(s)}').to_string()の戻り値はct02_Bonn_store_1: supply_Bonn_store_1 <= open_Bonnのような文字列なので:の後のsupply_Bonn_store_1 <= open_Bonn制約式だけを抜き出すために使っています。
スラック値はslack_valuesで取得します。

わかりやすくするためにDataFrameにしています。

制約2のスラック表示

import pandas as pd
def _extract_after_colon(s): return s[s.find(':') + 1:] if ':' in s else ''


warehouse = [w.id for w in warehouses for s in stores]
store = [s.id for w in warehouses for s in stores]
ct2n = [_extract_after_colon(
    warehouse_model.get_constraint_by_name(f'ct02_{str(w)}_{str(s)}').to_string())
    for w in warehouses for s in stores]
ct2s = [warehouse_model.get_constraint_by_name(f'ct02_{str(w)}_{str(s)}').slack_value
     for w in warehouses for s in stores]

pd.DataFrame({"warehouse": warehouse,
              "store": store,
              'ct': ct2n,
              'slack': ct2s})
              'slack': ct2s})

開設した倉庫のフラグとSupplyのフラグの差です。開設した倉庫には１が立ち、Supplyも1になったら、その差であるスラックは0になります。倉庫Bonnと店舗4のスラックは供給されたので０になります。BonnとBordeaux、Brusselsまでを表示しています。

制約3のスラック表示

同様に制約３のスラックをDataFrameで表示しています。

制約3のスラック表示

warehouse = [w.id for w in warehouses]
ct3n = [_extract_after_colon(
    warehouse_model.get_constraint_by_name(f'ct03_{str(w)}').to_string())
    for w in warehouses]
ct3s = [warehouse_model.get_constraint_by_name(f'ct03_{str(w)}').slack_value
     for w in warehouses]

pd.DataFrame({"warehouse": warehouse,
              'ct': ct3n,
              'slack': ct3s})

ParisとRomeは許容量をすべては満たしていないので1になっています。

結果のグラフ化

開設コストと供給コストの割合

開設コストと供給コストの割合を円グラフにします。
決定式の開設コストopening_cost_objと供給コストsupply_cost_objを[opening_cost_obj, supply_cost_obj]にしてpieチャートを書いています。

開設コストと供給コストの割合

if env.has_matplotlib:
    import matplotlib.pyplot as plt
    %matplotlib inline

def display_costs(costs, labels, colors):
    if env.has_matplotlib:
        plt.axis("equal")
        plt.pie(costs, labels=labels, colors=colors, autopct="%1.1f%%")
        plt.show()
    else:
        print("warning: no display")

display_costs(costs=[opening_cost_obj, supply_cost_obj],
              labels=["Opening Costs", "Supply Costs"],
              colors=["gold", "lightBlue"])

倉庫ごとの供給コストの割合

倉庫ごとの供給コストの割合を円グラフにします。
if supply_vars[w,s].solution_value >= 0.9で、倉庫と店舗の供給の有無supply_varsにフラグが立っている倉庫と店舗を見つけて、sum(supply_costs[w][s.id - 1] で輸送コストの合計を出しています。
なお、supply_varsの決定変数はBinaryですが、誤差のために１にならないことがあるので、solution_value >= 0.9でフラグが立っているかをチェックしています。

倉庫ごとの供給コストの割合

# 倉庫ごとの供給コストの割合
supply_cost_by_warehouse = [sum(supply_costs[w][s.id - 1]  
                                for s in stores if supply_vars[w,s].solution_value >= 0.9)
                             for w in warehouses]
wh_labels = [w.id for w in warehouses]
display_costs(supply_cost_by_warehouse, wh_labels, None)

パレートフロンティアの探索

ここまでは開設コストと供給コストの合計を最小化しましたが、ここでは開設コストと供給コストのそれぞれの最小値を確認してみます。

開設コストを優先して最小化する

開設コストを最小化してから、供給コストを最小化します。
set_multi_objectiveを使います。
senseは最小化minか最大化maxを決めます。
exprsには決定式のリストを入れます。ここではtotal_opening_cost, total_supply_costを入れています。
prioritiesで、total_opening_costとtotal_supply_costのどちらが優先されるのかを決めます。大きい値の方が優先されます。ここではtotal_opening_costが1、total_supply_costが0なので、total_opening_costが優先されることになります。
そして再度warehouse_model.solve()で求解しています。

開設コストを優先して最小化する

warehouse_model.set_multi_objective(
    sense="min", exprs=[total_opening_cost, total_supply_cost], priorities=[1, 0])
ok = warehouse_model.solve()
opening1 = total_opening_cost.solution_value
supply1 = total_supply_cost.solution_value

print("total opening cost=%g" % opening1)
print("total supply cost=%g" % supply1)
print("toltal cost is %g" % (opening1 + supply1))

results.append((opening1, supply1))

開設コストの絶対最小値は 120 であり、この値に対して最小供給コストは 352 であることがわかります。最初の試行では、合計を最小化して (140,293) を取得しましたが、開設コストについては120 が達成できる最良の値であることがわかりました。供給コストは 352 になり、合計コストは 472 になります。これは、合計を最適化したときに見つかった 433 という値よりも確かに大きいです。

結果

total opening cost=120
total supply cost=352
toltal cost is 472

供給コストを優先して最小化する

次は供給コストを最小化してから、開設コストを最小化します。
priorities=[0, 1]となりました。total_opening_costが0、total_supply_costが1なので、total_supply_costが優先されることになります。

供給コストを優先して最小化する

warehouse_model.set_multi_objective(
    sense="min", exprs=[total_opening_cost, total_supply_cost], priorities=[0, 1])
ok = warehouse_model.solve()
opening2 = total_opening_cost.solution_value
supply2 = total_supply_cost.solution_value

print("total opening cost=%g" % opening2)
print("total supply cost=%g" % supply2)
print("toltal cost is %g" % (opening2 + supply2))

results.append((opening2, supply2))

ここでは、(152, 288) が得られ、供給コストは 352 から 288 に下がりますが、開設コストは 120 から 152 に上がります。合計は 288+152 = 440 であり、合計の最適値 433 よりも大きいです。

結果

total opening cost=152
total supply cost=288
toltal cost is 440

パレート図

合計の最適化結果(Sum)、開始コスト優先(Opening_first)、供給コスト優先(Supply_first)の3 つの点を (開始コスト、供給コスト) 散布図にします。

パレート図

print (results)
def display_pareto(res):
    if env.has_matplotlib:
        plt.cla()
        plt.xlabel('Opening Costs')
        plt.ylabel('Supply Costs')
        colors = ['g', 'r', 'b']
        markers = ['o', '<', '>']
        nb_res = len(res)
        pts = []
        for i, r in enumerate(res):
            opening, supply = r
            p = plt.scatter(opening, supply, color=colors[i%3], s=50+10*i, marker=markers[i%3])
            pts.append(p)
        plt.legend(pts,
               ('Sum', 'Opening_first', 'Supply_first'),
               scatterpoints=1,
               loc='best',
               ncol=3,
               fontsize=8)
        plt.show()
    else:
        print("Warning: no display")

display_pareto(results)

完成pythonコード

参考

The Warehouse Problem

CPLEXサンプル warehouseを読み解く #OPL - Qiita

https://qiita.com/spssfun2017/items/c29c7ac9a06493f556d4
同じ問題のOPL言語版です。

Warehouse location - IBM Documentation

CPLEXサンプルを読み解く記事一覧