ppt

輸送問題に対する
主双対内点法
小崎敏寛
2011 7/21 22 RIMS2011
1
目的：輸送問題を数値的に解く
→主双対内点法の適用
•
•
•
•
•
線形計画問題に対する主双対内点法
輸送問題の定式化
輸送問題に対する内点法
内点法の解から整数解を得る手法
結論と今後の課題
2011 7/21 22 RIMS2011
2
1-1標準形の線形計画問題
主問題(P)
min cT x
s.t. A x  b
x  0 変数 xはｎ次元
双対問題(D)
max bT y
s.t. AT y  z  c 変数 yは m次元
z0 Ax  b
最適条件
変数 zはｎ次元
AT y  z  c
x0
z0
Xz  0
2011 7/21 22 RIMS2011
3
1-2 フィージブル内点法
Ax  b
AT y  z  c
x0
z0
センターパス
みたす
強くみたす
wk
wk+1
Xz   e
w  ( x, y , z )
最適解
各μにつき近似解を得て，μ↓0とすれば，
等式不等式系が最適条件に近づくので，最適解が得られるはず
主問題と双対問題を同時に考えるので主双対内点法
理論で使われる
2011 7/21 22 RIMS2011
4
1-3 インフィージブル内点法1
w  ( x, y , z )
w :ニュートン方向
Ax  b  Ax
wk
AT y  z  c  AT y  z
Zx  Xz   e  Xz
xT z
ただし   [0,1],  
n
w k 1  w k   cw
wk+1
ただし  c  0.99995,   [0,1]
最適解
実行可能領域の外から最適解に近づく
数値的に最適解を得られる
2011 7/21 22 RIMS2011
5
1-4 インフィージブル内点法2
• 計算時間の大半は次の方程式を解くのにか
かる
AXZ-1ATΔy=AXZ-1(c-ATy-z)-γμAZ-1e+b
Δz=-ATΔy+c-ATy-z
Δx=-XZ-1Δz+γμZ-1e-x
eは要素が全て1のベクトル
X:=diag(x) Z:=diag(z)
2011 7/21 22 RIMS2011
6
1-5 インフィージブル内点法3
• ステップサイズα+
min{ 1, max{ : x  x  0, z  z  0}}
• 非負であるようにステップサイズをきめる
2011 7/21 22 RIMS2011
7
2-1 輸送問題1
min
a1
供給
シンク
a2
b2
工場
a3
cij
M=3
N=2
ij
xij
(i , j )
b1
ソース
c
s.t. x11  x12
 a1
x21  x22
需要
 a2
x31  x32  a3
消費地
x11
 x21
x12
 x31
 x22
 b1
 x32  b2
xij  0
• 仮定：  ai   b j
i
j
• 仮定：全てのiとjの間にアークが存在する
• 制約が一つ冗長なので一番下の制約を取り除く
2011 7/21 22 RIMS2011
8
2-2 輸送問題2
min
c
ij
xij
(i , j )
s.t. x11  x12
 a1
x21  x22
 a2
x31  x32  a3
x11
 x21
 x31
 b1
xij  0
• MN変数M+N-1制約の
線形計画問題
• 実行可能解
c11
x11
c21
x21
c31
x31
b1
c12
x12 a1
c22
x22 a2
c32
x32 a3
b2
xij  ai b j /(  ai )
が存在する
i
2011 7/21 22 RIMS2011
9
2-3 輸送問題と他の問題の関係
（割当問題）（輸送問題）（多品種輸送問題）
輸送問題で，M=Nとして，ai , bj =1として，
変数を0-1に限定．
多品種輸送問題で品種を一つにしたもの．
2011 7/21 22 RIMS2011
10
2-4 輸送問題に対する内点法
の既存研究
• 変数に対する非負制約だけでなく，上限制約
のある定式化
←実行可能解が存在しないことがある
• ニュートン方向の計算に共役勾配法を使用
a1=3
u11=1
u12=1
2011 7/21 22 RIMS2011
11
3-1 輸送問題のニュートン方向
ij  xij / zij
• AXZ-1AT
 : (ij )

を の一番右の列を
M
取り除いた行列
＝
N-1
コレスキー分解
O((M+N-1)3)
M
計算量はO(MN)
N-1
密だとすると O(m2n)=O((M+N-1)2MN)
2011 7/21 22 RIMS2011
12
3-2 ステップサイズの決定（LP）
ステップサイズの決め方に問題がある
正値性をみたしていれば，
ニュートン方向には
できるだけ進んだ方がよいと考える
Infeasibilityはステップサイズが1以下
なら大きいほど改善する
min{ 1, max{  : x  x  0}
x2
Δx1
x
α
座標ごとに計算する
  i  min{ 1, max{  i : xi   i xi  0}
2011 7/21 22 RIMS2011
Δx
Δx2
x1
13
3-3 ステップサイズと残差の関係
Zx  Xz  e  Xsより
xTz
( x  x)T ( z  z )
γxTz
n
 xT z (1   (1   ))   2  xi zi
i 1
考えない
A( x  x)  bより
A( x  x)  Ax   ( Ax  b)
α
Ax
b
α
2011 7/21 22 RIMS2011
14
3-4MATLABでの実装で注意したこと
• diagを使用しない
→sparse(1:M*N,1:M*N,x./z)
• ステップサイズの計算の時に，先に領域を確
保しておく
2011 7/21 22 RIMS2011
15
3-5 数値実験結果
自作の内点法
linprog(lipsol)
反復数
時間(秒)
反復数
時間(秒) 時間(秒)
4
7.80×10-2
4
0.0160
0.0310
50
6
3.10×10-2
5
0.125
4.53
100
6
0.686
5
0.500
62.0
500
4
8.42
6
23.8
×
M=N
10
2011 7/21 22 RIMS2011
simplex
16
4 整数解の求解
• 以降では供給ベクトル，需要ベクトルは整数
ベクトルと仮定する．
• 内点法では，問題が退化している場合は，基
底解に収束するとはかぎらない．
→整数解が得られない
• 単体法では計算時間がかかる．
2011 7/21 22 RIMS2011
17
5-1 解の性質
• 1 センターパスの極限は強相補解
• 2 理論上の内点法で生成される点列は強相
補解に収束する[Güler and Ye1993]
• 3 実際の内点法で生成される点列は強相補
解に収束するのか？
2011 7/21 22 RIMS2011
18
5-2 解の性質2
• M=N=3の輸送問題
• x*+z*
x*z*
0.0129
0.0210
0.3315
0.3644
0.6556
0.6137
0.3315
0.3644
0.6667
0.3757
1.0018
0.4014
0.6556
0.6137
1.0018
0.4014
1.3426
0.3193
×10-9
強相補解が得られている
2011 7/21 22 RIMS2011
19
5-3 退化の種類
（主）退化（primal degeneracy）
主変数xの基底変数で値が0のものが存在
双対退化（dual degeneracy） [Gülerら1991]
双対変数zの非基底変数で0のものが存在
→対応する主変数xの非基底変数で正のもの
が存在
superbasic variable[Maros 2003]
（主変数xで）非基底変数で正の値をとるもの
2011 7/21 22 RIMS2011
20
5-4 内点法で得られた解
• 輸送問題（M=N=3）
x*=
0.0086
0.3321
0.6592
0.3321
0.6667
1.0012
0.6592
1.0012
1.3395
z*=
0.8630
0.5803
0.4941
0.5803
0.2975 ×10-6
0.2114
0.4941
0.2114
0.1253
2011 7/21 22 RIMS2011
x*=Θ(1)
z*=Θ(μ*)
21
6-1 整数解の求解
• 主双対内点法の解から基底解を得るアルゴリズム
(Megiddo 1988, Gülerら1991 )
添え字1をx>0の座標,添え字2をxとzが0の座標，添え
字3をz>0の座標とする．A=[A1 A2 A3], x=[x1 x2 x3],
c=[c1 c2 c3]とする．
• A1t=0をみたす非ゼロのtを求める．c1Tt=0がなりた
つ(近似的に)．
• x1の非負性を保つように，レシオテストをする．
• 0になったaｊをA1からA2に移す．
• A1の列ベクトルが独立になれば終了.A1xB=bを解き
xBを得れば整数解が得られる．
2011 7/21 22 RIMS2011
22
6-2 数値実験結果
• 輸送問題(M=N=9)
81変数の線形計画問題
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 3 0 0
0 0 0 0 0
0 0 0 4 0
0 2 0 0 0
1 0 0 0 5
0
0
0
0
0
0
0
6
0
0
0
0
0
0
2
3
0
2
0
2
0
0
2
4
0
0
0
2011 7/21 22 RIMS2011
1
0
3
4
0
0
0
0
1
行列形式の解x*
23
7 結論と今後の課題
• 輸送問題の係数行列の疎性を利用した計算
と（線形計画問題に対する）ステップサイズの
決定により効率的に最適解を得られた．
• より大きな輸送問題に対して，整数解を得る
アルゴリズムを考え，実装する．
2011 7/21 22 RIMS2011
24
参考文献
最近のサーベイ
• L. T. Guardia and G. B. Lima, Interior Point Methods for Linear
Multicommodity Flow Problems, 2nd International Conference on
Engineering Optimization, 2010.
内点法の収束
• O. Güler and Y. Ye, Convergence behavior of interior-point algorithms,
Mathematical Programming, Volume 60, Numbers 1-3, 215-228, 1993.
基底解の求解
• N. Megiddo, On finding primal- and dual-optimal bases, Research Report
RJ 6328, IBM Almaden Research Center, San Jose, California, 1988.
• O. Güler, D. den Hertog, C. Roos, T. Terlaky and T.Tsuchiya, Degeneracy in
interior point methods for linear programming, Report 91-102, 1991.
Superbasic variable
• I. Maros, Computational Techniques of the Simplex Method, Kluwer
Academic Pub, 2003.
MATLABの本
• 上坂吉則，MATLABプログラミング入門，牧野書店，2006.
2011 7/21 22 RIMS2011
25
• Δy（の近似）を求めるのに，直接法でなく反復
法を使う．
• SOR法(Successive Over-Relaxation，逐次
加速緩和法)を適用するが，3倍時間がか
かったので諦めた．
2011 7/21 22 RIMS2011
26
ステップサイズの決定（提案）
(SOCP,SDP)
(１)変数が二次錐の直積で表されている場合
  pi  min{ 1, max{ i : xi  i xi  0}
i
(2)SDPで変数がブロック対角の場合
各ブロックごとにステップサイズを決定する


pi
 min{ 1, max{  i : xi   i xi 0}
2011 7/21 22 RIMS2011
27

Download Report