ロジスティックネットワーク設計システム MEta Logistic network Optimization System

簡易システムビデオ

MELOS (MEta Logistic Optimization System)

MEta Logistic network Optimization System-GF (Green FIeld)

print("flow=", x.X )
print("transfer price", t.X)
for c in Country:
    print(f"Profit of {c}", Profit[c].X)
flow= 20000.0
transfer price 12.5
Profit of A 89999.99999999999
Profit of B 0.0
model.Params.NonConvex =2 #非凸の二次も解けるようする!
model.optimize()
Changed value of parameter NonConvex to 2
   Prev: -1  Min: -1  Max: 2  Default: -1
Gurobi Optimizer version 9.1.0 build v9.1.0rc0 (mac64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 4 rows, 13 columns and 11 nonzeros
Model fingerprint: 0x8543791d
Model has 3 quadratic constraints
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  QMatrix range    [1e+00, 1e+00]
  QLMatrix range   [1e+00, 2e+01]
  Objective range  [5e-01, 1e+00]
  Bounds range     [2e+04, 2e+04]
  RHS range        [0e+00, 0e+00]
  QRHS range       [2e+04, 1e+05]

Continuous model is non-convex -- solving as a MIP.

Presolve removed 1 rows and 1 columns
Presolve time: 0.00s
Presolved: 20 rows, 15 columns, 46 nonzeros
Presolved model has 4 bilinear constraint(s)
Variable types: 15 continuous, 0 integer (0 binary)

Root relaxation: unbounded, 5 iterations, 0.00 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0  postponed    0               -          -      -     -    0s
     0     0  postponed    0               -          -      -     -    0s
     0     2  postponed    0               -          -      -     -    0s
*  166   474              34    -89291.66667          -      -   1.0    0s
*  213   474              34    63242.130004          -      -   0.9    0s
*  308   474              33    108194.50543          -      -   0.7    0s
*  404   474              33    123928.57136          -      -   0.6    0s
* 1010   867              36    126709.65465          -      -   2.4    0s
* 1012   867              37    126727.74194          -      -   2.4    0s

Explored 11130 nodes (45579 simplex iterations) in 0.20 seconds
Thread count was 16 (of 16 available processors)

Solution count 6: 126728 126710 123929 ... -89291.7

Model is unbounded
Best objective 1.267277419356e+05, best bound -, gap -
for (i,j) in x:
    print("flow=",i,j, x[i,j].X )
    print("transfer price",i,j, t[i,j].X)
for c in Country:
    print(f"Profit of {c}", Profit[c].X)
flow= A B 18750.0
transfer price A B 7.516129032258065
flow= A C 20000.0
transfer price A C 7.516129032258065
Profit of A 0.0
Profit of B 97161.29032288078
Profit of C 111638.7096774193
126728/113600
1.1155633802816902

はじめに

ここでは,ロジスティクス・ネットワーク設計問題に対する包括的モデルを示す. モデルの目的は,単位期間(通常は年)ベースのデータをもとに, ロジスティクス・ネットワークの形状を決めることにある. モデルを求解することによって得られるのは,倉庫の設置の是非, 地点間別の各製品群の単位期間内の総輸送量, 工場別の各製品群の単位期間内の総生産量である.

ロジスティクス・ネットワーク設計モデルでは, 以下の意思決定項目に対する最適化を同時に行う.

  • 各製品群をどの工場でどれだけ生産するか?
  • 各製品群をどの倉庫(配送センターや中継拠点の総称)で保管するか?
  • 各顧客群の各製品群の需要を,どの地点から運ぶか?
  • 倉庫をどこに新設するか?(または移転,閉鎖するか?)
  • 複数の倉庫の候補地点からどれを選択するか?

このような意思決定は,単に現状のロジスティクス・ネットワーク の見直しを行う場合だけでなく,以下の状況においても有効に用いられる.

  • 吸収・合併の後のロジスティクス・ネットワークの再編成
  • 新製品投入時の意思決定(どこで製造して,どのような流通チャネルで流すか)
  • ロジスティクスにおける戦略的提携(ロジスティクス・パートナーシップ)

例題

ここで用いるサンプルデータは, データ生成モジュールのgenerate_cust, generate_prod, generate_demandなどの諸関数を用いて生成したものである.

以下,順に説明する.

顧客データ

品が最終的に消費される場所を顧客(customer)とよぶ. 顧客は,対象となるシステムに依存して,小売店であったり個々の家庭であったりする. モデルとして扱う場合には,同じ施設からサービスされると考えられる顧客の集まりを 1つの顧客に集約することができる. ここでは,個々の顧客を同じサプライ・チェインでサービス を行っても差し支えない程度に集約したものを顧客群(customer group)とよぶこともあるが,以下では単に顧客と呼ぶ.

通常は,需要の大きい顧客や特殊なサービスを要求する顧客は1つの顧客として扱い, 残りの顧客を地区ごとに集約して扱うことが多い. これは,同じ地区に属する顧客への配送は,(大口や特殊な顧客以外は) 同一のデポから行われる場合が多いことに起因する. 個々の顧客の数は数万から数十万におよぶこともあるが, 集約後の顧客群の数は,多くても数百になるようにする.

顧客群に集約するための簡便法として郵便番号を用いる方法もあるが, できれば1つの顧客群に含まれる顧客数および需要量が,なるべくバランスするように することが望ましい. ここでは,幾つかのクラスタリング手法 (Weizfeld法, $k$-means法, 階層的クラスタリング法, $k$-メディアン法) を用いる.

例題で用いる顧客データの列は, 名前 (name) と緯度・経度 (lat, lon; latitude, longitudeの略) である.

cust_df = pd.read_csv(folder + "Cust.csv", index_col=0)
cust_df.head()
name lat lon
id
1 札幌市 43.06417 141.34694
2 青森市 40.82444 140.74000
3 盛岡市 39.70361 141.15250
4 仙台市 38.26889 140.87194
5 秋田市 39.71861 140.10250

製品データ

ロジスティクス・ネットワーク内を流れるものを製品(product)とよぶ. 製品の種類は,対象とするロジスティクス・システムにもよるが,通常は膨大なものになるが, ストラテジックレベルの意思決定モデルを扱う場合には,通常同じ特性をもつ製品を集約して扱う. 製品を集約したものを,個々の製品と区別するために製品群(product group)とよぶこともある. 通常,製品群の数が数十程度になるように集約する. また,ABC分析をして,影響の少ない製品を除いて最適化を行うこともある.

製品データは,以下の列をもつが,ロジスティクス・ネットワーク設計モデルで使うのは,名称 (name)と重量 (weight)だけである.

  • name: 製品の名称
  • weight: 製品の重量 (単位はkg)
  • volume: 製品の容量 (単位は 𝑚3 )
  • cust_value: 顧客上での製品の価値
  • dc_value: 倉庫上での製品の価値
  • plnt_value: 工場における製品の価値
  • fixed_cost: 製品を生産する際の段取り(固定)費用
prod_df = pd.read_csv(folder + "Prod.csv", index_col=0)
prod_df.head()
name weight volume cust_value dc_value plnt_value fixed_cost
0 A 2 0 7 1 1 14
1 B 5 0 5 1 1 14
2 C 1 0 5 1 1 19
3 D 3 0 5 1 1 17
4 E 1 0 10 1 1 18

倉庫データ

最終需要地点および供給地点以外の場所で,かつ製品の製造・加工を行わない地点を総称して 倉庫(warehouse)とよぶことにする. 実際には,倉庫は流通センター(DC: Distribution Centerの略),デポ,クロスドッキング地点,港,空港,ハブなど様々な形態をとる.

倉庫の役割と機能には,以下のものが考えられる.

  • 調達(購入)の規模の経済性のため. 大量購入に伴う割引が行われるとき,大量に購入した製品を一時的に保管しておく際に,倉庫の保管機能が用いられる.

  • 混載による規模の経済性のため.複数の小売店(スーパーなど)で消費される製品を, 直接製造業者(メーカー)から輸送すると,小ロットの輸送のため費用が増大する. メーカーから車立て輸送を行うためには,複数の小売店の需要をあわせる必要があるが, そのための中継地点として倉庫の仕分け機能が用いられる. メーカーから運ばれた複数の製品は,倉庫で方面別に仕分けされて,小売店に運ばれる. この際,複数の製品を混載して輸送するので,やはり小ロット輸送による費用の増大を避けることができる. この仕分け機能に特化した倉庫をクロスドッキング地点(cross docking point)とよぶ. クロスドッキング地点は,保管機能をもたず,メーカーから小売店が発注した量と同じ量を受け取り, 即座に仕分けして小売店に輸送する.

  • リスク共同管理効果を得るため. 複数の小売店の在庫を,集約して倉庫に保管することによって,安全在庫量を減らすことができる.これをリスク共同管理(risk pooling)効果とよぶ. この場合には,倉庫は,リスク共同管理によってサプライ・チェイン全体での在庫を減らすために用いられる.

  • 顧客へのサービスレベル向上のため. たとえば,コンビニエンスストアの配送センターは,弁当が約$3$時間以内で配達可能なように設置される. これによって,$1$日に数回の弁当の配送が可能になり,顧客に出来たての弁当を届けることが可能になる. このように,倉庫は小売店(コンビニエンスストア)へのリード時間を減少させるために用いられる.

  • 単一ソース条件. 顧客(特に小売店)に対して,複数の倉庫からの配送を行うことを許さないという条件が付加されることがある. これは,荷捌き場の制限や,荷物の受け入れを複数回行う手間を省くためであるが,そのためには,複数のメーカーの製品を倉庫で積み替えて 運ぶことが必要になる.

  • 製品に付加価値をつけるため. 最終製品への分化を,なるべくサプライ・チェインの下流で行うことによって, 顧客サービスを向上させ,安全在庫を減少させることができる. これを遅延差別化(delayed differentiation, postponement)とよぶ. 遅延差別化の方法として, 工場で最終製品まで製造してまうのではなく,顧客に近い倉庫で,製造の最終工程を行うことがあげられる. たとえば,値札やラベルを付けたり,簡単な加工を行うことが代表例である.

倉庫データは,倉庫の配置可能地点を表し,以下の列をもつ.

  • name: 倉庫名称
  • lb: 容量下限
  • ub: 容量上限
  • fc: 固定費用
  • vc: 変動費用
  • lat: 緯度
  • lon: 経度
dc_df = pd.read_csv(folder + "DC.csv", index_col=0)
dc_df.head()
name lb ub fc vc lat lon
0 札幌市 0.0 5732551.52 10037 0.432510 43.06417 141.34694
1 青森市 0.0 5732551.52 10235 0.414573 40.82444 140.74000
2 盛岡市 0.0 5732551.52 10908 0.414802 39.70361 141.15250
3 仙台市 0.0 5732551.52 10072 0.136525 38.26889 140.87194
4 秋田市 0.0 5732551.52 10767 0.029622 39.71861 140.10250

工場データ

製品の製造・加工する地点を総称して工場(plant)とよぶ. 実際には,工場は生産工場,部品工場,流通加工を行う倉庫など様々な形態をとる.

工場における製品の生産量上限は,別途定義するので,工場は名称(name)と緯度・経度情報 (lat,lon) だけをもつ.

plnt_df = pd.read_csv(folder + "Plnt.csv", index_col=0)
plnt_df.head()
name lat lon
0 Odawara 35.284982 139.196133
1 Osaka 34.563101 135.415129
2 Chiba 35.543452 140.113908

需要データ

ロジスティクス・ネットワーク設計モデルでは,各顧客・製品に対する計画期間内の総需要量の情報が必要になる. 通常,計画期間は1年とし,年間の総需要を入力とする.これは,日々の需要を記録した需要データから計算される.

需要データは,以下の列をもつ.

  • date: 日付
  • cust: 顧客名
  • prod: 製品名
  • demand: 需要量
demand_df = pd.read_csv(folder + "demand.csv", index_col=0)
demand_df.head()
date cust prod demand sales
0 2019-01-01 札幌市 A 1 7
1 2019-01-01 札幌市 B 1 5
2 2019-01-01 札幌市 C 0 0
3 2019-01-01 札幌市 D 0 0
4 2019-01-01 札幌市 E 0 0

生産情報データ

計画期間内における工場における製品の生産量には上限がある.これを定義するのが生産情報データである. 使用する列は,工場名 (plnt),製品名 (prod),生産量上限 (ub)だけである.

plnt_prod_df = pd.read_csv(folder + "Plnt-Prod.csv", index_col=0)
plnt_prod_df.head()
plnt prod ub lead_time
0 Osaka G 5282 28
1 Osaka F 4649 29
2 Osaka J 12526 25
3 Osaka D 3313 26
4 Osaka E 17964 28

Weiszfeld法を実行する関数 weiszfeld_numpy と weiszfeld

施設数が1つの場合には、Weiszfeld法と呼ばれる反復法で近似解が得られることが知られている。 複数施設に対しては、k-means法と似た反復解法が考えられる。ここでは、重み付き複数施設配置問題に対するWeiszfeld法を構築する。

移動距離は大圏距離distanceを用いる。これはgeopyから読み込んであると仮定する。 例題においては,顧客の重みは、年間の製品需要量に製品の重量を乗じたものとする。

引数:

  • cust_df : 顧客データフレーム;緯度経度情報を利用する。
  • weight : 顧客の重み;重み付きの大圏距離の和を最小化するものとする。
  • num_of_facilities : 平面上に配置する施設の数を指定する。
  • epsilon : 誤差上限
  • max_iter : 反復回数の上限
  • seed: 乱数の種

返値:

  • X : 施設の緯度のリスト
  • Y : 施設の経度のリスト
  • partition : 顧客の割り当てられた施設の番号を格納したリスト
  • cost: 費用(重み付きの距離の総和)

weiszfeld_numpyはNumPyを用いた実装であり, weiszfeldは通常の実装である.

weiszfeld_numpy[source]

weiszfeld_numpy(cust_df, weight, num_of_facilities, epsilon=0.0001, max_iter=1000, seed=None)

Weiszfeld法; 複数施設の連続施設配置問題の近似解法 (NumPy version)

weiszfeld[source]

weiszfeld(cust_df, weight, num_of_facilities, epsilon=0.0001, max_iter=1000, seed=None)

Weiszfeld法; 複数施設の連続施設配置問題の近似解法

weiszfeld関数の使用例

cust_df = pd.read_csv(folder + "Cust.csv", index_col=0) #read customer data for using Lat/Lng
cust_df.head()
name lat lon
id
1 札幌市 43.06417 141.34694
2 青森市 40.82444 140.74000
3 盛岡市 39.70361 141.15250
4 仙台市 38.26889 140.87194
5 秋田市 39.71861 140.10250
prod_df = pd.read_csv(folder + "Prod.csv", index_col=0) #read product data (for using weight)
prod_df.head()
name weight volume cust_value dc_value plnt_value fixed_cost
0 A 2 0 7 1 1 14
1 B 5 0 5 1 1 14
2 C 1 0 5 1 1 19
3 D 3 0 5 1 1 17
4 E 1 0 10 1 1 18
total_demand = pd.read_csv(folder + "total_demand.csv")
total_demand.head()
Unnamed: 0 prod cust demand
0 0 A さいたま市 93.255495
1 1 A 京都市 165.453297
2 2 A 仙台市 62.170330
3 3 A 佐賀市 907.486264
4 4 A 前橋市 85.233516
weight_of_prod ={p:w for (p,w) in zip(prod_df.name, prod_df.weight) } #prepare weight dic.
weight_of_cust = defaultdict(int)
for row in total_demand.itertuples():
    weight_of_cust[row.cust] += weight_of_prod[row.prod]*row.demand 
#weight = [weight_of_cust[row.name] for row in cust_df.itertuples()] 
weight = [weight_of_cust[i] for i in weight_of_cust] 
# 関数呼び出し
#X, Y, partition, cost = weiszfeld(cust_df, weight, num_of_facilities = 10, epsilon=0.0001, max_iter = 10, seed=1)
#print(cost)

X, Y, partition, cost = weiszfeld_numpy(cust_df, weight, num_of_facilities = 10, epsilon=0.0001, max_iter = 10, seed=1)
print(cost)
11030261.63674016

反復Weiszfeld法を実行する関数 repeated_weiszfeld

Weiszfeld法を一定回数初期解を変えながら実行し,最も良い結果を返す関数

引数:

  • cust_df : 顧客データフレーム;緯度経度情報を利用する。
  • weight : 顧客の重み;重み付きの大圏距離の和を最小化するものとする。
  • num_of_facilities : 平面上に配置する施設の数を指定する。
  • epsilon : 誤差上限
  • max_iter : 反復回数の上限
  • numpy: NumPyを用いる場合にTrue
  • seed: 乱数の種

返値:

  • X : 施設の緯度のリスト
  • Y : 施設の経度のリスト
  • partition : 顧客の割り当てられた施設の番号を格納したリスト
  • cost: 費用(重み付きの距離の総和)

repeated_weiszfeld[source]

repeated_weiszfeld(cust_df, weight, num_of_facilities, epsilon=0.0001, max_iter=1, numpy=True, seed=None)

TODO: Multicore (threading)で高速化

%time
X, Y, partition, cost = repeated_weiszfeld(cust_df, weight, 
            num_of_facilities = 10, epsilon=0.0001, max_iter = 30, numpy=False, seed=2)
print(cost)
CPU times: user 2 µs, sys: 1e+03 ns, total: 3 µs
Wall time: 6.91 µs
9674288.59942924

連続施設配置の可視化関数

引数:

  • cust_df : 顧客データフレーム
  • X : 施設の緯度のリスト
  • Y : 施設の経度のリスト
  • partition : 顧客の割り当てられた施設の番号を格納したリスト
  • weight: 顧客の重み(この大きさにあわせて点を描画する.)

返値:

  • fig: Plotlyの図オブジェクト

show_optimized_continuous_network[source]

show_optimized_continuous_network(cust_df, X, Y, partition, weight=None)

連続施設配置の可視化関数

年間(計画期間)の総需要を計算する関数 make_total_demand

引数:

  • demand_df: 需要データフレーム
  • start: 計画期間の開始日
  • finish: 計画期間の終了日
  • num_of_days: 推定したい総需要の日数(既定値は365;すなわち年間需要量の推定)

返値:

  • 総需要を入れたデータフレーム

make_total_demand[source]

make_total_demand(demand_df, start='1900/01/01', finish='2050/12/31', num_of_days=365)

年間(計画期間)の総需要を計算する関数

demand_df = pd.read_csv(folder + "demand.csv")
total_demand_df, demand_cust_df, demand_prod_df  = make_total_demand(demand_df, start="2019/01/01", finish="2050/12/31", num_of_days=365)
weight = demand_cust_df.demand.values
#total_demand_df = make_total_demand(demand_df)
total_demand_df.head()
prod cust demand
0 A さいたま市 3650.000000
1 A 京都市 28635.909091
2 A 仙台市 2123.636364
3 A 佐賀市 136908.181818
4 A 前橋市 5906.363636

顧客の集約

顧客の集約はロジスティクス・ネットワーク設計モジュールのrepeated_weiszfeldでできる.scikit-learnにも幾つのかのクラスタリング法が実装されているが,顧客の重み(需要量)を考慮したものとしてKmeans法とその変形があげられる.

これらの手法を用いて適切な顧客の集約を行うことを考える.

kmeans法を実行する関数 kmeans

kmeasn法は顧客の重み付き2乗距離の和を最小化する. 本来の重み付き距離の総和ではないが,高速なので大規模問題例に簡易的に使うことができる.

引数:

  • cust_df : 顧客データフレーム;緯度経度情報を利用する。
  • weight : 顧客の重み;重み付きの大圏距離の和を最小化するものとする。
  • num_of_facilities : 平面上に配置する施設の数を指定する。

返値:

  • X : 施設の緯度のリスト
  • Y : 施設の経度のリスト
  • partition : 顧客の割り当てられた施設の番号を格納したリスト
  • cost: 費用(重み付きの距離の総和)

kmeans[source]

kmeans(cust_df, weight, num_of_facilities=1, batch=True)

kmeans関数の使用例

repeated_weiszfeldより悪いが,高速である.

%%time
X, Y, partition, cost = repeated_weiszfeld(cust_df, weight, num_of_facilities = 8, epsilon=0.0001, max_iter = 3, numpy=True, seed=1)
print(cost)
490929454.9189146
CPU times: user 973 ms, sys: 1.3 ms, total: 974 ms
Wall time: 973 ms

階層的クラスタリング関数 hierarchical_clusterning

hierarchical_clusterning[source]

hierarchical_clusterning(cust_df, weight, durations, num_of_facilities=2, linkage='average')

階層的クラスタリング関数 hierarchical_clusterning

durations,  distances, node_df  = compute_durations(cust_df)
X, Y, partition, cost = hierarchical_clusterning(cust_df, weight, durations, num_of_facilities = 10, linkage="complete")
print(cost)
50725609599.77609

k-メディアン問題を解くための関数 solve_k_median

引数:

  • cust_df : 顧客データフレーム;緯度経度情報を利用する。
  • weight : 顧客の需要量
  • cost: 移動費用(地図で求めたdurationsを入力する.)
  • num_of_facilities : 平面上に配置する施設の数を指定する。
  • max_iter =100: 最大反復数
  • max_lr =0.01: 最大学習率(ステップサイズ)
  • moms=(0.85,0.95): 慣性項の下限と上限
  • convergence=1e-5: 収束判定のためのパラメータ
  • lr_find = False: 学習率探索を行う場合 True
  • adam = False: Adam法によるパラメータの調整を行う場合 True,上界と下界の差をもとに決める場合 False

返値:

  • X : 施設の緯度のリスト
  • Y : 施設の経度のリスト
  • partition : 顧客の割り当てられた施設の番号を格納したリスト
  • best_ub: 費み付きの費用の総和
  • lb_list: 下界値の変化を入れたリスト
  • ub_list: 上界値の変化を入れたリスト
  • phi_list: 学習率探索の場合の学習率の変化を入れたリスト

solve_k_median[source]

solve_k_median(cust_df, weight, cost, num_of_facilities, max_iter=100, max_lr=0.01, moms=(0.85, 0.95), convergence=1e-05, lr_find=False, adam=False)

k-メディアン問題を解くための関数 solve_k_median

solve_k_medianの使用例

#weight=np.random.random(len(cust_df))*10
durations, distances, node_df  = compute_durations(cust_df)
X, Y, partition, best_ub, lb_list, ub_list, phi_list = solve_k_median(cust_df, weight, durations, num_of_facilities=50, 
                                 max_iter=10000, max_lr=50., moms=(0.85,0.95), convergence=1e-5,  lr_find = False, adam=True)
0: 10.000, 0.00000, 0.00000, 47.00

k-メディアンの学習率探索の結果を可視化する関数 plot_k_median_lr_find

引数:

  • phi_list: 学習率設定のパラメータ $\Phi$ の変化を入れたリスト
  • lb_list: 下界の変化を入れたリスト

返値:

  • Plotlyの図オブジェクト

plot_k_median_lr_find[source]

plot_k_median_lr_find(phi_list, lb_list)

k-メディアンの学習率探索の結果を可視化する関数 plot_k_median_lr_find

plot_k_median_lr_find関数の使用例

fig = plot_k_median_lr_find(phi_list, lb_list)
plotly.offline.plot(fig);

k-メディアン最適化の探索を可視化する関数 plot_k_median

引数:

  • lb_list: 下界の変化を入れたリスト
  • ub_list: 上界の変化を入れたリスト

返値:

  • fig: 下界と上界の変化を可視化した図オブジェクト

plot_k_median[source]

plot_k_median(lb_list, ub_list)

plot_k_median関数の使用例

fig = plot_k_median(lb_list, ub_list)
plotly.offline.plot(fig);

k-メディアン法で得られた結果