BTST 法の 3 次元への拡張について
M2014SS007 中山莉奈
指導教員：鈴木敦夫

1

はじめに本 研 究 で は ，2 次 元 の Big Triangle Small Trianglemethod (BTST 法) [3] を 3 次元に拡張した Big Tetrahedron Small Tetrahedron method (BTeSTe 法) を新たに提案する．

需要点数が多い場合 BSSS 法より計算時間が長くなることである．理由は，BTST 法は初めに需要点 n をドローネ分割することで作成した 2n − 5 の三角形の計算時間が必要だからである．BSSS 法を 3 次元に拡張した解法として

Big Cube Small Cube method (BCSC 法)[7] も提案されている．

2 次元の BTST 法は目的関数が非凸である配置問題の本研究の目的は，以下の 3 つの 3 次元の配置問題につ最適解を分枝限定法を用いて求める解法のことである．一

いて BTeSTe 法を適用し，その有効性を確かめることで般的に，目的関数が非凸である配置問題は局所最適解がある．

複数存在する．そのため，このような配置問題は大域的最• 迷惑施設配置問題 [5]

適解を見つけることが困難である．それに対して，BTST• WAR 問題 [6]法は大域的最適解を比較的簡単に見つけることができる．

Drezner and Suzuki[3] では迷惑施設配置問題と WAR 問
題 (Weber problem with Attraction and Repulsion) についてその有効性が示されていた．また，Drezner andBrimberg[2] では真円度問題についてその有効性が示されていた．しかし，従来の BTST 法では 2 次元の配置問題• 真球度問題
ここで，3 つの 3 次元の配置問題について簡単に記述する．
前提として，各配置問題は空間に分布した需要点とその重みを持つこととする．
迷惑施設配置問題とは，新しく施設を配置するとき危険しか解くことができない．そこで更に多くの問題が解けるや被害が最小となる地点を決定する配置問題のとこであよう，BTST 法を 3 次元に拡張した BTeSTe 法を提案する．例として，原子力発電所やゴミ処理場の配置問題が挙る．また，この解法を MATLAB 上で実現し，計算機実験げられる．

を行うことでその有効性を確認する．
ここで，BTeSTe 法の概略を記す．BTeSTe 法は 2 次元

WAR 問題とは，Weber 問題の重みに負の値が加わった配置問題のことである．また，Weber 問題は配置場所からの場合と同様に分枝限定法を用いる．制約条件として，最既存施設へのユークリッド距離の重み付きの総和が最小と適解は需要点の凸包内にあることとし，扱う問題は非凸なる地点を決定する配置問題のことである．

線形最適化問題である．BTeSTe 法の手順は，初めに制真球度問題とは，部品などの表面から幾つかの点を選び約領域内をドローネ分割を用いて需要点を重なりのないその座標をデータとしたとき，球とみなしてよいかを判定四面体に分割する．次に四面体ごとに目的関数の Lowerする問題のことで鋼球の製品検査などに用いられている．

Bounds(LB) を計算する．そして，計算した LB と実行可能解から求めた Upper Bound(U B) を比較する．このU B より値の大きい LB を持つ四面体を考慮対象から除く

本研究では，半径とユークリッド距離間の差の絶対値の総和を最小化することを考える．
本研究では 3 次元の BTeSTe 法を上記に挙げた 3 つのことで最適解の存在する領域を狭めていく．これにより， 配置問題に適用し，MATLAB 上で計算機実験を行うこと従来は厳密解を求めることが困難であった非凸線形最適化 でその有効性を確かめる．また，BCSC 法との比較を計算問題に対しても厳密解を求めることができる．

BTST 法より以前に，非凸線形最適化問題に対して厳密解を求めることができる解法が提案されていた．その解法

機実験を用いて行う．

2BTeSTe 法のアルゴリズム

は Big Square Small Square technique (BSSS 法) [4] で2 次元の BTST 法を 3 次元に拡張したときのアルゴリある．BSSS 法は BTST 法の三角形の代わりに四角形をズムを説明する．本研究は最小化問題を仮定する．各需用いて解を求める解法である．BSSS 法は初めに実行可能

要点の座標 Xi = (xi , yi , zi )(i = 1, ..., n) と重み wi (i =領域である四角形を 1 つ置き，その四角形を分割してい1, ..., n) は所与とする．そして制約条件として，最適解はくことで最適解の存在する領域を狭めていく．この解法と需要点の凸包内にあることとし，扱う問題は非凸線形最適

比べ BTST 法の利点は，解が凸包内に存在するかどうか化問題である．目的関数の Lower Bounds の計算方法はを調べる必要がないことである．理由はドローネ分割を用

Drezner[1] に詳しく記載されている．

いるからである．BSSS 法は初めに全ての需要点を含むよ

1. 需要点を用いた 3 次元ドローネ分割

うな四角形で区切るため，凸包外も解の検索範囲に含まれ需要点はドローネ分割を用いて重なりのない四面体にる．そのため，凸包内に解を限定する場合は四角形が凸包分割する．ドローネ分割する理由は，実行可能領域全内に存在するかどうかを調べる必要がある．逆に欠点は，体がカバーできるからである．ここで，ドローネ分割

図 3 八面体の分割方法

図 1 凸多角形のドローネ分割の一例

れる各四面体は正四面体に近い形に分割できる．の一例を図 1 に示す．

そのため，極端に細長い形の四面体が減り計算時2. Upper Bound (U B) の計算間の短縮につながると考える．このようにして 1各四面体の重心を用いて評価を行う．そして，それらつの四面体を 8 つの小さな四面体に細分割する．の中で最小の目的関数値を暫定解とし，U B と設定

その後，細分割した 8 つの小さな四面体の重心でするの目的関数値と LB をそれぞれ求める．目的関数

3. 各四面体の Lower Bounds (LB) の計算

値は U B より値が小さくなったとき U B の値を4. U B/(1 + ϵ) と LB の比較 (ϵ：十分小さい値)その目的関数値に更新する．LB は 8 つの小さなのとき四面体のうち 1 つをリスト先頭の分割前の四面その LB を持つ四面体を考慮対象から除外する体と入れ替えダウンヒープを行う．残り 7 つの小（a）LB > U B/(1 + ϵ)（b）LB < U B/(1 + ϵ)

さな四面体は 1 つずつリストの最後に追加しアッのときその四面体内に解が存在する可能性があるため，プヒープを行うことでリストを管理する．LBmin

リストにその四面体の情報を格納する．ここでいの値が U B/(1 + ϵ) より大きくなるまでその比較うリストとは，考慮対象である四面体の情報を管を繰り返す．

理するものである．格納する情報は四面体の頂点番号と LB の値である．リストはヒープを用い配置問題の定式化と計算機実験の結果3て管理し，常に一番上の情報がリスト内で最小のこ こ で は 3 つ の 配 置 問 題 に つ い て 説 明 し ，こ れ ら

LB （LBmin ）を持つ四面体になるようにする．
LBmin の値が U B/(1 + ϵ) より大きいとき，アルゴリズムを終了する．このとき，最適解は U B となる．LBmin の値が U B/(1 + ϵ) より小さいとき，その四面体を 8 つの小さな四面体に細分割すの 配 置 問 題 に BTeSTe 法 ，BCSC 法 を 適 用 し 計 算 機
実 験 を 行 な っ た 結 果 を 記 す ．使 用 す る 計 算 機 は Win-をとり，元の頂点 1 つと 3 つの中点から成る四面dows7 Intel(R)Core(TM) 2 Duo CPU P8700 @2.53GHz
2.53GHz，ソフトウェアは MATLAB R2013a である．また，計算機実験は ϵ = 10−6 で行い，需要点の座標はランダムに与える．そして，10 個の問題を作成し計算した結果の平均を記す．2 つの解法は制約条件として，最適解は需体を 4 つ作成する (図 2 の I，II，III，IV)．

要点の凸包内にあることと仮定する．本研究は最小化問題

る．四面体の細分割方法は，はじめに各辺の中点を取り扱う．
次に，表記方法を定義しておく．{X1 , ..., Xn } は需要点n の集合で，それらの需要点 i の座標は Xi = (xi , yi , zi ) と示し，配置場所 X の座標は X = (x, y, z) と示す．このとき配置場所 X と需要点 i 間のユークリッド距離を di (X)とする．図形 s の頂点 Tj の座標は XsTj で示す．Xsc は図形 s の重心座標を示す．j は図形の頂点番号を示し，BTeSTe 法のとき j = 1, ..., 4，BCSC 法のとき j = 1, ..., 8となる．
図 2 頂点と中点からなる四面体
次に，残りの中点のみで構成された八面体を 4 つ

3.1

迷惑施設配置問題

うち，1 番短い対角線を挿入することで 4 つの四

迷惑施設配置問題は各需要点の迷惑度の大きさを wi ≥0(i = 1, ..., n) とする．この迷惑度は距離の 2 乗に反比例面体に分割する (図 3 の V，VI，VII，VIII)．ここして減少すると仮定する．このとき，迷惑度の総和が最小で 1 番短い対角線を用いる理由は，計算時間の短となる点が解となる．そのため，目的関数は式 (1) のよう縮を計るためである．１番短い対角線から構成さになる [3]．の四面体に分割する．分割方法は 3 本の対角線のF (X) =n∑

wi2 (X)
di=1 i(1)

迷惑施設配置問題は配置場所が需要点と等しい (X =Xi ) とき，ユークリッド距離が di (X) = 0 となる．そのた

このとき F + (X)，F − (X) はどちらも凸になる．WAR 問
題の目的関数は F (X) = F + (X) − F − (X) の最小化となる．しかし，この 2 つの凸関数間の差は一般的に局所最適解を含む．そのため，WAR 問題の目的関数は凸になるとは限らない．

め，迷惑度が ∞ となり，目的関数の一部が ∞ になる．したがって，迷惑施設配置問題の目的関数は非凸である．

3.1.1

3.2.1

WAR 問題の LB の計算方法は Drezner and Suzuki[3]Lower Bounds の計算方法

迷惑施設配置問題の LB の計算方法は Drezner andSuzuki[3] の考え方を用いる．d2i (X) = x，d2i (Xsc ) = a とし，関数 y = 1/x の (a, 1/a) における接線の方程式を利用する．つまり，式 (2) のように下から値をおさえる．

11
1x
1≥ − 2 (x − a) = (2 − )xa aaaの考え方を用いる．式 (6) のように F + (X) の接平面を用いて下から値をおさえる．

F + (XsTj ) ≥ F (XsTj )∑ { xsc − xi
= F + (Xsc ) +wi[xsTj − xsc ]di (XsTj )+i∈I}ysc − yizsc − zi+[ysTj − ysc ] +[zsTj − zsc ]di (XsTj )di (XsTj )s(6)(2)

式 (2) を利用して図形 s の重心での目的関数値を式 (3) のように下からおさえる．

よって，図形 s の重心での目的関数値は式 (7) のように下F (Xsc ) ≥ min LsTj

(3)j

からおさえることができる．

()
di (XsTj )wi
LsTj =2−d2 (Xsc )d2i (Xsc )i=1 in∑F (Xsc ) = F + (Xsc ) − F − (Xsc ) ≥ min LsTj
jj

このとき，各図形 s の LB は下記のようになる．

(s = 1, ..., S)

LBs = min LsTj ,(s = 1, ..., S)

j

3.1.2

計算機実験の結果

3.2.2

MATLAB を用いて算法を実現し，計算した迷惑施設配置問題の結果を表 1 に記す．迷惑度はランダム (wi ≤ 1)に与えた．

点数

102050100
3.2

反復回数

BTeSTe2210.62179.1
1994.45441.7

MATLAB を用いて算法を実現し，計算した WAR 問題の結果を表 2 に記す．重みはランダム (−1 ≤ wi ≤ 1) に
表 2 WAR 問題の計算結果

計算時間 (秒)BCSC
3969.57710.68783.87803.1

BTeSTe
3.73
5.05
8.44
15.35

計算機実験の結果与えた．

表 1 迷惑施設配置問題の計算結果
需要(7)LsTj = F (XsTj ) − F − (XsTj )このとき，各図形 s の LB は下記のようになる．

LBs = min LsTj ,Lower Bounds の計算方法

需要BCSC
219.70733.591499.321485.72

点数102050100

反復回数

BTeSTe327.6
708.5
1401.02368.9

計算時間 (秒)

BCSC
114.3
197.0
283.5
480.5

BTeSTe
0.65
1.27
2.52
4.54

BCSC
4.93
15.39
28.4287.89

WAR 問題
真球度問題

WAR 問題は Weber 問題の各需要点の重み wi が正負の値をとる配置問題のことである．Weber 問題の目的関数

3.3

は既存施設へのユークリッド距離の重み付きの総和の最小のことである．本研究は式 (8) のように半径 r とユーク化のため，目的関数は式 (4) のようになる [3]．

リッド距離 di (X) の差の絶対値の総和が最小となる点をF (X) =n∑求める．

wi di (X)

(4)

i=1

F (X) =

式 (4) は式 (5) のように重みの正負で分けることができる．
需要点が正の重みを持つときを I + ，負の重みを持つときを I − と仮定する．

F + (X) =∑i∈I +wi di (X),

真球度問題は球と見なしてよいかどうかを判定する問題

F − (X) =∑n∑|di (X) − r|

(8)

i=1

しかし，この問題は di (X) − r の計算より正負の値が混在するため解くことが困難である．そこで，di (X) のメディ{−wi di (X)}i∈I −(5)アン距離を半径 r と仮定し，Drezner and Brimberg[2] で提案された順序メディアン問題を用いることで真球度問題を式 (9) のように定式化する．

n∑

F (X) =|di (X) − median {di (X)|i ∈ n}|

表 4 球面内に予め図形を入れた時の計算結果(9)

需要i=1点数a(q) (X) は di (X) の 距 離 を 短 い 順 に 並 び 替 え た も の(a(1) (X) ≤ ... ≤ a(n) (X)) とする．このとき，式 (9) は式(10) のように置き換えることができる．
n∑
F (X) =|di (X) − median{di (X)}|i=1=∑ (a(q) (X) − median{di (X)}q>n/2

∑ (−=2a(q) (X) −∑

3.4

BCSC
3.29
2.67
1.80
2.19

考察が有効であることがわかった．また，解の大まかな場所がdi (X)(10)

分かる問題は，その部分を含むような図形を初めに入れ，そこから計算を行った方が計算時間が早くなることがわ

この 2 つの凸関数間の差は一般的に局所最適解を含む．そのため，真球度問題の目的関数は凸になるとは限らない．

3.3.1

計算時間 (秒)BTeSTe
3.95
4.28
2.56
2.58

BTeSTe 法，BCSC 法を用いた計算機実験の結果より，解が凸包上になる可能性が高い配置問題は BTeSTe 法の方

i∈nq>n/2

BCSC
782.0
567.5
275.1237.3))
a(q) (X) − median{di (X)}

q≤n/2

∑

10
20
50
100

反復回数

BTeSTe1706.71639.3714.1
517.8

Lower Bounds の計算方法

かった．その場合，BTeSTe 法と BCSC 法では大きな差はなかった．

4おわりに

真 球 度 問 題 の LB の 計 算 方 法 は Drezner and

本研究は BTST 法を 3 次元に拡張した BTeSTe 法を新
Brimberg[2] の 考 え 方 を 用 い る ．式 (11) の よ う にたに提案した．そして，この解法を3 つの 3 次元の配置問∑
2 q>n/2 a(q) (X) の接平面を用いて下から値をおさえる．
題に適用し，計算機実験を行うことで実用的な時間内で問(x(q) , y(q) , z(q) ) は a(q) (Xsc ) に対応する需要点 (q) の座標
題が解けることを確認した．

を示す．

U (XsTj ) = 2{∑(

a(q) (Xsc ) +(xsc − x(q) )(xsTj − xsc )a(q) (Xsc )参考文献[1] Z. Drezner, A general global optimization approachfor solving location problems in the plane, Journal(ysc − y(q) )(ysTj − ysc ) (zsc − z(q) )(zsTj − zsc ) )}++of Global Optimization, 37, 305–319, 2007.
a(q) (Xsc )a(q) (Xsc )(11) [2] Z. Drezner, J. Brimberg, Fitting concentric circles tomeasurements, Mathematical Methods of Operations
したがって，図形 s の重心での目的関数値は式 (12) のよResearch, 79, 119–133, 2014.
うに下からおさえることができる．
n{}∑[3] Z. Drezner, A. Suzuki, The big triangle small trianF (Xsc ) ≥ min U (XsTj ) −di (XsTj )(12)jgle method for the solution of non-convex facility loi=1
cation problems, Operations Research, 52, 128–135,
このとき，各図形 s の LB は下記のようになる．
2004.
n{}
∑LBs = min U (XsTj ) −di (XsTj ) , (s = 1, ..., S) [4] P. Hansen, D. Peeters, D. Richard, J.F. Thisse,q> n2

j3.3.2i=1

The minisum and minimax location problems revisited,Operations Research, 33, 1251–1265, 1985.

計算機実験の結果

MATLAB を用いて算法を実現し，計算した真球度問題の結果を表 3 に記す．半径 (r = 20) から p1/5（但し，p は0 から 1 の一様分布に従う乱数）を引いた結果の座標データを用いている．
表 3 真球度問題の計算結果
需要点数

反復回数

計算時間 (秒)102050100

1970.51200.7894.6789.8

4.61
3.20
3.41
4.23

次に球面内に予め四面体または立方体を入れ，その図形から解法を用いて計算を行った結果を表 4 に記す．

[5] 中山莉奈, 鈴木敦夫, BTST 法の 3 次元への拡張について, 日本 OR 学会 2015 年春季研究発表会, 1-D-4,2015.
[6] R. Nakayama, A. Suzuki, Z. Drezner, T. Drezner,
The Big Tetrahedron Small Tetrahedron Method forThree Dimensional Location Problems, INFORMS.
Philadelphia, 2015-11.
[7] A. Schobel, D. Scholz, The big cube small cube solution method for multidimensional facility location
problems ,Computers & Operations Research, 37,
115–122, 2010.