ＲでＯＲ：動的計画法

ナップサック問題
生産在庫問題
輸送問題
スケジューリング問題
関数の最小値問題
巡回セールスマン問題

動的計画法

動的計画法は、数理モデルの解を数値的に求める際に、多段階決定過程に置き換えて、効率的に解を計算する方法である。 Smith, "Dynamic Programming" (1994) の演習問題をRを使って解く。動的計画法の漸化式は合成積のように、「たすき掛け」の計算が基本なので、ちょっとしたプログラミングが必要になる。プログラミングの手順は

動的計画法ＤＰの漸化式を書いて
ステージも状態と同じように添え字を使って表し（たとえば、f_n(S) => f[n,S]）
漸化式を数式で表せばよい。

C言語を使っても同じことだけれど、Rの場合は一々コンパイルしなくて良いので、デバッグがスムースに行く。Cのようなプログラミング言語に慣れていると、添字の動かし方や、条件判定の仕方が異なるので、注意が必要である。たとえば、

添え字は1以上しか使えないので、漸化式に出てくるゼロをそのまま生かすためには、i(k) = k+1 という関数を用意して、添字 k の代わりに i(k) とするとよい。
「for(i in 1:n)」の「1:n」は「１から n までの等差数列」という意味で、n>0 ならば公差は１、n<=0 の場合は公差は -1 となることに注意する。
複合条件文は個別の条件式をすべてチェックするので、たとえば、「i=0; if(i>0 & a[i]>0」は a[0] が定義されていないのでエラーになる

ナップザック問題

重さと価値が与えられているいくつかの品物と、重量制限のある容器が与えられているとき、品物の中から適当に選んで容器に入れて、その価値を最大にしたい。品物の組合せと、価値の最大値はいくつか、という問題をナップザック問題という。ナップザックになるべくたくさんものを詰め込む、というイメージから名付けられたものであろう。

ナップザック問題（重複なし）

問題　次のような4種類の品物が与えられている。重量の上限が与えられたとき、価値を最大にする品物の組合せを求めよ。ただし、品物はすべて１個ずつしかないものとする。

品物の番号 k	1	2	3	4
重量 w(k)	2	3	4	6
価値 v(k)	8	12	15	19

品物の数が k 種類で、重量制約が n の場合の最適解（価値の合計）を f(k,n) とすると、次のような漸化式が成り立つ。

f(k,n) = max{ f(k-1, n-w(k))+v(k), f(k-1, n) } f(k,n) = f(k-1, n)　　　　　　　　　　　　　　　　n < w(k) の場合はこちら

この式をそのままプログラムにすれば良い。このとき、2番目の添字は０まで動くので、そのままではエラーになるため、添字を1ずらす関数を定義しておくと、デバッグが楽になる。

> i <- function(x) {x+1}			# 添字を一つずらす関数

データベクトルを作成し、m を品物の数、T を重量上限とする。

> w <- c(2,3,4,6)				# 「重さ」ベクトル
> v <- c(8,12,15,19)				# 「価値」ベクトル
> T <- sum(w)					# 最大の重さ
> m <- length(v)				# 品物の個数

最初に、(m+1) x (T+1) の大きさの作業用テーブルを用意して、すべての要素に０を代入する。

> F <- matrix(0, m+1, T+1)			# 作業用テーブル　(m+1) x (T+1)

品物が１種類しかない場合、上限値が w(1) 以上ならば、最適解は v(1) である。

> F[1,i(w[1]):i(T)] <- v[1]			# 品物が一つしかない場合

k 種類ある（k > 1）場合は、漸化式にしたがって、k 番目の品物を含めた場合と含めない場合のどちらか大きいものを選べば良い。
max() 関数をベクトルの各要素に適用するためには、いったん行列にして apply() 関数を使うと、for 文を使わずに、ベクトル計算が出来る。

> for(k in 2:m) {
+   F[k,i(0):i(w[k]-1)] = F[k-1,i(0):i(w[k]-1)]
+   F[k,i(w[k]):i(T)] = apply(rbind(F[k-1,i(w[k]):i(T)],F[k-1,i(0):i(T-w[k])]+v[k]), 2, max)
+ }

ここで、rbind() 関数は行ベクトルを行列にまとめる関数、apply(A,2,max) は、行列の各列（2番目の引数の"2"）について、最大値（3番目の"max"）を適用（計算）する関数である。その結果、次が得られる。

> colnames(F) = 0:T
> F
     0 1 2  3  4  5  6  7  8  9 10 11 12 13 14 15
[1,] 0 0 8  8  8  8  8  8  8  8  8  8  8  8  8  8
[2,] 0 0 8 12 12 20 20 20 20 20 20 20 20 20 20 20
[3,] 0 0 8 12 15 20 23 27 27 35 35 35 35 35 35 35
[4,] 0 0 8 12 15 20 23 27 27 35 35 39 42 46 46 54

結果の行列をそのまま表示させると、列番号が１から割り振られる。列は価値に対応していて、左端は０と表示して欲しいので、colnames 関数を使い、列に名前を付与して表示してある。

ある重量制限のもとで最適値を達成する品物の組合せは、漸化式で作った数表 f を見れば求めることが出来る。
実際、f で列に同じ数字があるのは、その容量（列の値）ならば、品物が増えても、増えた品物を使わなくても最大価値が達成されていること、そうでないならば、新たに追加された品物を取り込むことによって、価値の大きい組合せが達成させること、を意味する。
この問題設定で、容量11の問題の最適解は、次のようにして求められる。f(4,11) ≠ f(3,11) なので、最初にアイテム４を取り込むと、残りの容量は5になる。次に f(3,5) = f(2,5) なので、アイテム３は無視する。同様にして、f(2,5)≠f(1,5) なのでアイテム２を取り込み、最後は f(1,2)≠0なので、アイテム１も取り込む。結局、アイテム1,2,4を取り込んで、合計価値は39となる解が得られた。

これらの考察から、品物が j 個で上限が k の場合の組合せ g(j,k)(j=1,...,m) を求めるには次のようにすればよい。f(j,k) と f(j-1,k) が等しければ、品物 j は使わないので、f(j-1,k) を達成する組合せと同じ、等しくない場合は、品物 j が最適解を達成するために必要なので、f(j-1,k-w(j)) を達成する組合せに品物 j を加えたもので最適解を達成する。

f(m,k) を達成する品物の組合せ（j を使うならば１、さもなければ０とする m 次元のベクトル）を g(1:m,k) とすれば、g を計算するアルゴリズムは次で与えられる。

g(1:m, 0:T) = 0 とする。k=1 から T まで次を実行する。ka = k とする。
j = m とする。
もし、f(j,ka) = 0 あるいは f(j,ka) = f(j-1,ka) ならばステップ6へ。
g(j,k) = 1 とする。ka - w(j) をあらたな ka とする。
ka = 0 ならば、ステップ7へ。
j-1 をあらたな j とし、j>1 ならばステップ3へ。
ka > 0 で f(1,ka) > 0 ならば、g(1,k) = 1

R はベクトル計算が得意なので、ベクトル演算を念頭に置いたアルゴリズムに書き換えると、例えば次のようになる。

f(m,n] ≠ f(m-1,n) ならば1、さもなければ0(n=0,...,T)というベクトルを g(m,) とする。
h = g(m,) * w(m) とする。
k=m-1,...,2 について４、５を実行する。
f(k,n-h(n)) ≠ f(k-1,n-h(n)) ならば1、さもなければ0(n=0,...,T)というベクトルを g(k,) とする。
h = h + g(k,) * w(k) とする。
f(1,n-h(n)) ≠ 0 ならば1、さもなければ0(n=0,...,T)というベクトルを g(1,) とする。

以上の説明に基づいて、各品物の重量ベクトル w と価値ベクトル v を与えて、品物の個数と重量上限値の様々な組合せに対して最大価値を計算し、その最大価値を生み出す品物の組合せを求める関数 knapsackproblem() を作る。

knapsackproblem = function(w, v) {
  i = function(x) {x+1}
  T = sum(w)
  m = length(v)	
  # 最適値行列
  F = matrix(0, m, T+1)
  F[1,i(w[1]):i(T)] = v[1]
  for(k in 2:m) {
    F[k,i(0):i(w[k]-1)] = F[k-1,i(0):i(w[k]-1)]
    F[k,i(w[k]):i(T)] = apply(rbind(F[k-1,i(w[k]):i(T)],F[k-1,i(0):i(T-w[k])]+v[k]), 2, max)
  }

  # アイテムの組み合わせ
  G = matrix(0, m, T+1)
  G[m,] = as.numeric(F[m,] != F[m-1,])
  h = G[m,] * w[m]
  for(k in (m-1):2) {
    G[k,] = as.numeric(F[k,i(0:T-h)] != F[k-1,i(0:T-h)])
    h = h + G[k,]*w[k]
  }
  G[1,] = as.numeric(F[1,i(0:T-h)] != 0)
  colnames(F) = colnames(G) = 0:T
  return(list(value=F, item=G))
}

# 実行例
> w <- c(2,3,4,6)				# 「重さ」ベクトル
> v <- c(8,12,15,19)				# 「価値」ベクトル
> knapsackproblem(w, v)
$value
     0 1 2  3  4  5  6  7  8  9 10 11 12 13 14 15
[1,] 0 0 8  8  8  8  8  8  8  8  8  8  8  8  8  8
[2,] 0 0 8 12 12 20 20 20 20 20 20 20 20 20 20 20
[3,] 0 0 8 12 15 20 23 27 27 35 35 35 35 35 35 35
[4,] 0 0 8 12 15 20 23 27 27 35 35 39 42 46 46 54

$item
     0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
[1,] 0 0 1 0 0 1 1 0 0 1  1  1  1  0  0  1
[2,] 0 0 0 1 0 1 0 1 1 1  1  1  0  1  1  1
[3,] 0 0 0 0 1 0 1 1 1 1  1  0  1  1  1  1
[4,] 0 0 0 0 0 0 0 0 0 0  0  1  1  1  1  1

最適解は $value 行列の右下と、$item 行列の右端列を見れば良い。この場合は、全て取り込んで 54 の価値を得る。この $value 行列、$item 行列は容量１５以下の全ての問題に対する最適解も与えている。例えば、容量１１の問題は１１の列を見れば良い。それによると、アイテム1,2,4を取り込んで価値 39のナップザックを作ることができる。

ナップザック問題（重複あり）

（42ページの問題）

重量ベクトルと価値ベクトルが与えられているとき、
容量Ｔのナップザックに最大価値を詰め込むにはどうしたらよいか。
同じものを複数選んでも良い

プログラム例

w = c(6,4,3,2)		# 重量
v = c(19,15,12,8)	# 価値
f = matrix(0, 5, 21)
for(n in 1:4) {
  for(T in 0:20) {              
    f[n+1,T+1] = f[n-1+1,T+1]
    if(T/w[n] < 1) next          ### for 文はCと違い、1回は必ず実行される
    for (j in 1:(T/w[n])) {
       f[n+1,T+1] = max(f[n+1,T+1], j*v[n] + f[n-1+1,T-j*w[n]+1])
    }
  }
}
f

ナップザック問題（複数制約問題）

（47ページの問題）

重量ベクトルと容量ベクトルと価値行列が与えられている。価値は非線形
容量Ｔのナップザックに重量／容量制約の下で最大価値を詰め込む
同じものを複数選んでも良い。

プログラム例

w = c(10,7,5,4,2)	# 重量
c = c(18,20,16,15,5)	# 容量
v = cbind(c(30,50,65,73,73,73,73), c(20,38,50,55,60,60,60),
          c(18,28,32,36,36,36,36), c(14,16,18,20,22,24,26),
          c(10,18,25,30,30,30,30))	# 価値（列ベクトルの結合）
f = array(0,dim=c(101,41,6))		# 3次元行列の定義
for(n in 1:5) {
	for(T in 0:40) {
		for(S in 0:100) {
			f[S+1,T+1,n+1] = f[S+1,T+1,n-1+1]
			if(min(S/c[n], T/w[n], 7) < 1) next   ### 価値行列が7個分まで指定
			for(i in 1:min(S/c[n], T/w[n], 7)) 
				f[S+1,T+1,n+1] = max(f[S+1,T+1,n+1], v[i,n] + f[S-i*c[n]+1,T-i*w[n]+1,n-1+1])
		}
	}
}
for(T in 40:35)
cat(f[101:96,T+1,6], "\n")

ページトップにもどる　索引にもどる

生産在庫問題

基本問題

（70ページの問題１）

需要予測が与えられている
保管費が与えられている（10までは250、それ以上は200）
生産する場合はセットアップコストがかかる
費用最小の生産計画を求めたい

プログラム例

D = c(11,15,19,24,28,11)                   ## 予測需要量（逆順）
iv = c(0, 250*(1:10), 2500+200*(1:sum(D))) ## 保管費
su = 3000                                  ## セットアップコスト

f = matrix(0,7,sum(D)+1)  ## 費用行列

Imax = 0
for(n in 1:6) {
	Imax = Imax + D[n]
	for(I in 0:Imax) { ### 計算が必要な領域
		if(I >= D[n]) {
			f[n+1,I+1] = iv[I-D[n]+1] + f[n-1+1,I-D[n]+1] ## 生産しなくても良い場合 
			if(Imax == I) next
			for(x in 1:(Imax-I)) {
				f[n+1,I+1] = min(f[n+1,I+1], iv[I+x-D[n]+1] + su + f[n-1+1,I+x-D[n]+1])
			}
		} else { 
			f[n+1,I+1] = f[n-1+1,1] + su
			if(D[n]-I+1 > Imax-I) next
			for(x in (D[n]-I+1):(Imax-I)) {
				f[n+1,I+1] = min(f[n+1,I+1], iv[I+x-D[n]+1] + su + f[n-1+1,I+x-D[n]+1])
			}
		}
	}
}
f

すべての可能な初期在庫量についてすべて計算しているので、メモリがかなり必要。少し考えれば、需要量の一部だけを繰り越しでまかない、不足分を生産する、という中途半端な方策は最適ではないことが分かるので、計算すべき初期在庫量は限られる（次のプログラム）。検算用に残しておく。

生産在庫問題（メモリ縮小版）

（70ページの問題１）

需要予測が与えられている
保管費が与えられている（10までは250、それ以上は200）
生産する場合はセットアップコストがかかる
費用最小の生産計画を求めたい

プログラム例

D = c(11,15,19,24,28,11)                   ## 予測需要量（逆順）
iv = c(0, 250*(1:10), 2500+200*(1:sum(D))) ## 保管費
su = 3000                                  ## セットアップコスト

f = matrix(0,7,7)  ## 費用行列
g = matrix(0,7,7)  ## 計算すべき初期在庫量
h = matrix(0,7,7)  ## 生産量（gで換算）

for(n in 1:6)      ## 計算に必要な各期の初期在庫量
	for(i in 1:n)
		g[n+1,i+1] = g[n+1,i-1+1] + D[n-i+1]

for(n in 1:6) {
	f[n+1,0+1] = su + iv[sum(D)+1]  ## 繰り越し在庫無し（必ず生産）
	h[n+1,0+1] = 1
	for(i in 0:(n-1)) {
		if(su + iv[g[n-1+1,i+1]+1] + f[n-1+1,i+1] < f[n+1,0+1]) {
			h[n+1,0+1] = i+1
			f[n+1,0+1] = su + iv[g[n-1+1,i+1]+1] + f[n-1+1,i+1]
		}
	}
	for (j in 1:n) {                ## 繰り越し在庫あり
		f[n+1,j+1] = iv[g[n-1+1,j-1+1]+1] + f[n-1+1,j-1+1]
		h[n+1,j+1] = 0
		for(i in (j+1):(n-1)) {
			if(n > j+1)
				if(su + iv[g[n-1,i]+1] + f[n-1+1,i+1] < f[n+1,j+1]) {
					f[n+1,j+1] = su + iv[g[n-1,i]+1] + f[n-1+1,i+1]
					h[n+1,j+1] = i-j
				}
		}
	}
}
f
h

生産在庫問題、その２

（70ページの問題4）

需要予測が与えられている
製造単価が与えられている（逓減、非線形）
1台以上は作らなければいけない、10台までしか作れない
保管費が与えられている（線形）
費用最小の生産計画を求めたい
必ず1台は生産しなければいけないので、繰り越し量に???

I_1 <= D_1-1,I_2 <= D_2+D_1-2,...

という制約が付く。したがって、生産量に関しても

I_{n}+x_{n}-D_{n} <= D_{n-1}+...+D_1-(n-1),x_{n} >= 1

という制約付き。

プログラム例

iv = 15
f = matrix(0,7,sum(D)) ## 評価関数
h = matrix(0,6,sum(D)) ## 最適政策

Imax = 0
for(n in 1:6) {
	Imax = Imax + D[n] - 1           ## 生産量、繰り越し量の上限値
	for(I in 0:Imax) {
		x0 = max(1, D[n]-I)            ## 最小限作る場合
		f[n+1,I+1] = x0*q[x0] + iv*(I+x0-D[n]) + f[n-1+1,I+x0-D[n]+1]
		opt = x0
		h[n,I+1] = opt
		if(x0 >= Imax-1-I) next        ## それ以上は作れない場合
		for (x in (x0+1):min(10,Imax-1-I)) {              ## 最小値の計算
			if(x*q[x] + iv*(I+x-D[n]) + f[n-1+1,I+x-D[n]+1] < f[n+1,I+1]) {
				f[n+1,I+1] = x*q[x] + iv*(I+x-D[n]) + f[n-1+1,I+x-D[n]+1]
				opt = x                    ## 最小値の更新
			}
		}
		h[n,I+1] = opt
	}
}
f
h

生産在庫問題、その３

（70ページの問題5）

毎月の需要予測が与えられている
月産生産能力とコストが与えられている（「残業」処理をするので、コストが折れ線）
保管費が与えられている（線形）
費用最小の生産計画を求めたい
残業代が高いときは前期に作っておいて繰り越した方が有利???

プログラム例

cost = function(n, I, x) {
	rr[n] * min(x,regular[n]) + ro[n] * max(0,x-regular[n]) + iv*(I+x-D[n])
}
widgets = function(D) {
	f = matrix(0,6,sum(D)+1) ## 評価関数
	h = matrix(0,5,sum(D)+1) ## 最適政策
	Imax = 0
	for(n in 1:5) {
		Imax = Imax + D[n]                           ## 生産量、繰り越し量の上限値
		for(I in 0:Imax) {
			x0 = max(0, D[n] - I)
			f[n+1,I+1] = cost(n,I,x0) + f[n-1+1,I+x0-D[n]+1]
			opt = x0
			h[n,I+1] = opt
			x9 = min(regular[n]+overtime[n], Imax-I)
			if(x0 >= x9) next                          ## それ以上は作れない場合
			for (x in (x0+1):x9) {                     ## 最小値の計算
				if(cost(n,I,x) + f[n-1+1,I+x-D[n]+1] < f[n+1,I+1]) {
					f[n+1,I+1] = cost(n,I,x) + f[n-1+1,I+x-D[n]+1]
					opt = x                                ## 最小値の更新
				}
			}
			h[n,I+1] = opt
		}
	}
	z = rep(f[6,1],6)
	z[5] = h[5,1]
	y = z[5] - D[5]
	for(i in 4:1) {
		z[i] = h[i,y+1]; y = y + z[i] - D[i]
	}
	return(rev(z))
} 

rr = c(16,18,18,15,15)
ro = c(22,24,25,20,20)
regulara = c(400,300,250,350,350)
overtimea= c(120,150,150,100,150)
iv = 3
#Da = c(300,400,350,400,300)
Da = c(300,400,350,200,300)
D = Da/10; regular = regulara/10; overtime = overtimea/10 

widgets(D)

ページトップにもどる　索引にもどる

輸送問題

階層化された規則的なネットワークの輸送問題

（21ページ例題）

都市間の距離行列が与えられている
出発地点から目標地点まで最短距離の経路を求めたい
都市が階層的に規則的に配置され、各パスの枝の本数は等しい

p = rbind(c(99,99,99,99,99,99,99,99,99,99),
          c(29,99,99,99,99,99,99,99,99,99),
          c(23,99,99,99,99,99,99,99,99,99),
          c(99,13,15,99,99,99,99,99,99,99),
          c(99,11,18,99,99,99,99,99,99,99),
          c(99,18,16,99,99,99,99,99,99,99),
          c(99,99,99,22,19,99,99,99,99,99),
          c(99,99,99,22,16,22,99,99,99,99),
          c(99,99,99,99,15,16,99,99,99,99),
          c(99,99,99,99,99,99,15,19,22,99)) ## 距離行列
r = cbind(c(1,1,1),c(2,3,3),c(4,5,6),c(7,8,9),c(10,10,10)) ## 先行関係
                        ## 行列構造を使うため、要素数を揃える必要がある
f = matrix(0,10,5)

for(n in 2:5) {
	for (i in r[,n]) {
		f[n,i] = p[i,i] + f[n-1,i]
		for(j in r[,n-1]) 
			f[n,i] = min(f[n,i], p[i,j] + f[n-1,j])
	}
}
f

輸送問題（一般のネットワーク）

（26ページ例題）

都市間の距離行列が与えられている
出発地点から目標地点まで最短距離の経路を求めたい
一般の道路網を対象としている

p = rbind(c( 0,99,99,99,99,99,99,99,99,99),
          c(29,99,99,99,99,99,99,99,99,99),
          c(23,99,99,99,99,99,99,99,99,99),
          c(99,13,15,99,99,99,99,99,99,99),
          c(99,11,18,99,99,99,99,99,99,99),
          c(99,18,16,99,99,99,99,99,99,99),
          c(99,99,99,22,19,99,99,99,99,99),
          c(99,99,99,22,16,22,99,99,99,99),
          c(99,99,99,99,15,16,99,99,99,99),
          c(99,99,99,99,99,99,15,19,22,99)) ## 距離行列
N = 10; f = matrix(0,N,N)
f[1,] = p[1,]; f[1,1] = 0 
				## f[n,i] : iから出発して、nステップで到着する最短距離
for(n in 2:N) { 
	for(i in 1:N) {
		f[n,i] = p[i,1] + f[n-1,1]
			for(j in 2:N) f[n,i] = min(f[n,i], p[i,j]+f[n-1,j])
	}
}
f

洗浄液と生産効率の問題

71ページ問題６

日産8個のラインＡでは１リットルの洗浄液が必要、しかしその内の0.5リットルは回収して翌日にも使える。
日産5個のラインＢでは１リットルの洗浄液が必要、しかしその内の0.75リットルは回収して翌日にも使える。
使える洗浄液に限りがある場合、どのラインを稼働させるのが良いか。

a = 8; b = 5; as = 2; bs = 1
f = matrix(0,6,20)
h = matrix(0,6,20)

for(n in 1:5) {
	for(x in 2:(n*as)) {
		if(a + f[n-1+1,x-as+1] > b + f[n-1+1,x-bs+1]) {
			h[n,x] = 1
			f[n+1,x+1] = a + f[n-1+1,x-as+1]
		} else {
			h[n,x] = 2
			f[n+1,x+1] = b + f[n-1+1,x-bs+1]
		}
	}
}
f
h

ページトップにもどる　索引にもどる

スケジューリング問題

処理時間と納期が与えられたn個のジョブを、総納期遅れ最小化でスケジュール

nCm = function(n,m) {
	c = 1;
	if(m == 0 || m == n) return(c)
	for(i in 1:m) c = c * (n-i+1) / i
}
nCmp = function(n,m,aa) {  ## パターンベクトルの生成
	if(m == 0) {
		cc = c(rep(0,n),aa); ##cat(cc,"\n")
		return(cc)
	}
	cc = c(rep(1,m)); if(m < n) cc = c(cc,rep(0,n-m)) 
	cc = c(cc ,aa); ##cat(cc,"\n")
	for(k in m:(n-1)) {
		bb = c(1); if(k < n) bb = c(bb,rep(0,n-1-k))
		cc = c(cc, nCmp(k,m-1,c(bb,aa)))
	}
}

d = c(108,174,140,153,45)
p = c(75,81,51,34,23)
#p=c(24, 8,24,31,18)
#d=c(53,29,39,50,69)
n = 5
base = c(); for(i in 1:n) base = c(base, 2^(i-1))
f = rep(0, 2^n)
for(m in (n-1):0) {
	z = matrix(nCmp(n,m,c()),,n,byrow=T) ## n個からm個を選ぶ
	for(i in 1:nrow(z)) {
		ID = sum(z[i,] * base)
		cst = 10000
		for(j in 1:n) {
			if(z[i,j] == 1) next
			m = sum(replace(z[i,],j,1) * base)
			cst = min(cst, max(sum(z[i,]*p)+p[j]-d[j], 0) + f[m+1])
		}
		f[ID+1] = cst
		cat(ID, cst, "\n")
	}
}

for(m in (n-1):0) {
	z = matrix(nCmp(n,m,c()),,n,byrow=T)
	for(k in 1:nCm(n,m)) {
		for(j in 1:5) {
			if(z[k,j] == 1) cat(j," ") else cat(0," ")
		}
		cat(":",f[sum(z[k,] * base)+1],"\n")
	}
}

ページトップにもどる　索引にもどる

関数の最小値問題

1変数関数の最小値を動的計画法で探索する。解の存在範囲を絞り込む方法

## フィボナッチ探索による関数の最小化問題の解法 ################
#    f は関数定義
#    [a, b] は初期点
#    err は収束したと見なす x の範囲

fibSearch = function(f,a,b,err) {
	phi = (1+sqrt(5))/2
	x = (b-a)/phi/phi + a
	y = (b-a)/phi + a
	c(a,x,y,b," ==> ",f(a),f(x),f(y),f(b))
	while(abs(x-y) > err) {
		if(f(x) < f(y)) {
			b = y; y = x; x = (b-a)/phi/phi + a
			cat(a,x,y,b," ==> ",f(a),f(x),f(y),f(b),"\n")
		} else {
			a = x; x = y; y = (b-a)/phi + a
			cat(a,x,y,b," ==> ",f(a),f(x),f(y),f(b),"\n")
		}
	}
	return((x+y)/2)
} 

f = function(x) {sin(x)+exp(x)-x}
a = -1; b = -0.79; ymin = f(-0.924)*0.999; ymax = f(-1)
(fibSearch(f,a,b,0.001))

注釈

関数の引数は以下の通り
- 「f」は関数名、最後から3行目に定義されている
- 「a,b」は区間の端点
- 「err」は誤差の許容範囲（最大誤差）
「phi」は黄金比、最初の2点（x<y）は黄金比を利用して決める
あとは、区間を[a,x]とするのか[y,b]とするのか、区間の幅が誤差の許容範囲に収まるまで繰り返す

グラフィカルプレゼンテーション付き

## フィボナッチ探索による関数の最小化問題の解法 ################
#    求解過程を視覚化するために、計算点をプロットする
#    f は関数定義
#    [a, b] は初期点
#    [ymin, ymax] はプロットするために用いる関数値の範囲
#    err は収束したと見なす x の範囲

fibSearch = function(f,a,b,ymin,ymax,err) {
	phi = (1+sqrt(5))/2
	A = a; B = b
	x = (b-a)/phi/phi + a
	y = (b-a)/phi + a
	plot(c(a,x,y,b),c(f(a),f(x),f(y),f(b)),xlim=c(A,B),ylim=c(ymin,ymax))
	while(abs(x-y) > err) {
		if(f(x) < f(y)) {
			b = y; y = x; x = (b-a)/phi/phi + a
			points(c(x,y),c(f(x),f(y)))
		} else {
			a = x; x = y; y = (b-a)/phi + a
			points(c(x,y),c(f(x),f(y)))
		}
	}
	return((x+y)/2)
} 

f = function(x) {sin(x)+exp(x)-x}
a = -1; b = -0.79; ymin = f(-0.924)*0.999; ymax = f(-1)
(fibSearch(f,a,b,ymin,ymax,0.01))

ページトップにもどる

巡回セールスマン問題

n 個の都市からなるネットワークで、どの都市の間も直接行き来ができるとしたとき、全ての都市を１回ずつ訪問してもとの都市の戻る最短経路を見つけよ、という問題。都市の数を n、都市 i と k の間の距離を d(i,k) とする。都市を 1,2,...,n と番号付け、 1,...,n の順列を J=(j₁,...,j_n) とすると(j₁ = j_n+1 = 1)、Σ_id(j_i, j_i+1) を最小化するような順列 J を見つけよ、と定式化できる。

すでに n-m 個の都市を回った後、いま都市 k にいるとしよう。今までに訪問した都市の集合を S として、f_m(S,k) を残りの都市を巡回するために要する最短時間とすると、k の次に訪問する都市までの時間とそこから先の残り時間の和を最小化すれば良いので、次の式が成り立つ。

f_m(S,k) = min_not(i∈S) { d(k,i) + f_m-1(S∪{i},i) }

初期値を f₀{1,...,n},k) = d(k,1) としてこの漸化式を計算することにより、この問題の最適解が f_n-1({1},1) で求められる。

プログラムは、上の漸化式を忠実になぞっただけ、ただし、f の下付き文字を付ける代わりに、再帰を使って実現している。

tsp = function(d) {
  tspsub = function(S,k) {
    if(length(S) == n) return(d[k,1])		## 全部回り終わったら最後は1へ戻る
    more = 1000000
    for(i in setdiff(2:n, S)) {
      temp = d[k,i] + tspsub(c(S,i),i)		## ｋからiへ移動
      if(temp < more) {
        more = temp
        index = i
      }
    }
    if(more+dist(S) < opt) {			## 最適解の更新（外部変数）
      opt <<- more + dist(S)
      SS <<- c(S, index)
    }
    return(more)
  }

  n = dim(d)[1]
  SS = numeric(n)
  opt = 1000000
  return(list(optimalValue=tspsub(1,1), optimalRoute=SS))
}

プログラムの注釈：length 関数はベクトルの要素数を返す。この if 文が再帰呼び出しの終了条件。
setdiff 関数は差集合を返す。したがって、for 文は、まだ S に含まれない都市全てについて所要時間を計算するために使われる。
それまでに見つかった最適解を更新した場合は、外部変数に経路を記録している。 <<- は永久代入演算子と言い、外部変数の値を変更する場合に使われる。

実行例を載せる。

> d = matrix(c(0,10,20,22,28,38,52, 10,0,12,31,35,40,32, 20,12,0,50,45,34,31, 22,31,50,0,14,26,55,
+              28,35,45,14,0,24,60, 38,40,34,26,24,0,25, 52,32,31,55,60,25,0), 7)
> tsp(d)
$optimalValue
[1] 138

$optimalRoute
[1] 1 2 3 7 6 5 4

> d = rbind(c( 0,47,56,48,33,43),
+           c(47, 0,43,30,20,44),
+           c(56,43, 0,24,53,53),
+           c(48,30,24, 0,46,30),
+           c(33,20,53,35, 0,47),
+           c(43,44,53,30,47, 0))	## 都市間距離行列

> tsp(d)
$optimalValue
[1] 193

$optimalRoute
[1] 1 5 2 3 4 6

ページトップにもどる