強化学習： ノート10

読んでいる本（出典）： 強化学習 : Richard S.Sutton, Andrew G.Barto, 三上 貞芳, 皆川 雅章 : 本 : Amazon.co.jp

前回：ノート9 ／ 次回：ノート11
目次：強化学習

今日読んだページ： 98～118ページ
以下、自分の解釈。

• 有限 MDP であるような強化学習問題の解法の1つである動的計画法（DP）は、最も基本的には以下。
  • 方策反復手法：
    • 反復方策評価で初期方策の価値関数を（近似的に）求める。
    • 求まった価値関数について greedy な方策、つまり、常に を選択するような新しい方策に更新する（これは元の方策よりよいことが保証される）。
    • これを繰り返すとどんどん方策をよくしていくことができる。
  • 価値反復手法：
    • 初期方策すら与えず、任意に初期化した価値関数を、個々の状態に対して現時点の価値関数で最適な手を選んで（ ＝ Bellman 最適方程式で）更新していく。
    • これを繰り返すとどんどん価値関数が最適価値関数に近づくので、じゅうぶん更新幅が小さくなったら求まった価値関数について greedy な方策を解とする。

---

2016-04-11 修正
103～105ページの、2つのレンタカー営業所でレンタカーを融通し合う例題をやってみる。
方策反復手法で実装して、0初期化した状態価値関数と方策から開始する。

• 「方策評価」→「方策改善」を1ループ実行した状態。

• 「方策評価」→「方策改善」を4ループ実行した状態 ＝ 最適方策。教科書105ページの図 と同じはず。この例題では 5 ループ目で方策の変更が起こらず、4ループ目で最適方策に達していることがわかる。

スクリプトは以下。

# 各営業所への要求台数と返却台数の確率分布 (在庫なしや駐車台数オーバー未考慮)
p.n.car.demand.raw <- list(office1=dpois(0:20, lambda=3), office2=dpois(0:20, lambda=4))
p.n.car.return.raw <- list(office1=dpois(0:20, lambda=3), office2=dpois(0:20, lambda=2))
# 配列の末尾に、21台以上要求/返却される確率を足し上げて付加しておく
for (office in 1:2) {
  p.n.car.demand.raw[[office]][[22]] <- 1.0 - sum(p.n.car.demand.raw[[office]][1:21])
  p.n.car.return.raw[[office]][[22]] <- 1.0 - sum(p.n.car.return.raw[[office]][1:21])
}

# 補助関数： 朝の時点で何台にしておくと夜何台になるか & いくら報酬がもらえるか
GetNextState <- function(office, n.car.morning) {
  p.n.car <- rep(0.0, 21) # 夜時点で0～20台になっている確率
  reward  <- rep(0.0, 21) # 上記のときの報酬 (＝ 10ドル * ちゃんと貸せた要求台数)
  
  # ありうる要求台数 (在庫なし考慮)
  n.car.demand   <- seq(0, n.car.morning, by=1)
  p.n.car.demand <- p.n.car.demand.raw[[office]][n.car.demand + 1]
  p.excess <- sum(p.n.car.demand.raw[[office]][(n.car.morning + 1 + 1):22]) # 在庫なしになる確率
  p.n.car.demand[[n.car.morning + 1]] <-  p.n.car.demand[[n.car.morning + 1]] + p.excess
  for (i.d in 1:length(n.car.demand)) {
    # ありうる返却台数 (駐車台数オーバー考慮)
    n.car.return.max <- min(20, max(0, 20 - n.car.morning + n.car.demand[[i.d]]))
    n.car.return   <- seq(0, n.car.return.max, by=1)
    p.n.car.return <- p.n.car.return.raw[[office]][n.car.return + 1]
    p.excess <- sum(p.n.car.return.raw[[office]][(n.car.return.max + 1 + 1):22]) # 駐車台数オーバーになる確率
    p.n.car.return[[n.car.return.max + 1]] <-  p.n.car.return[[n.car.return.max + 1]] + p.excess
    for (i.r in 1:length(n.car.return)) {
      n.car <- n.car.morning - n.car.demand[[i.d]] + n.car.return[[i.r]]
      p.n.car[[n.car + 1]] <- p.n.car[[n.car + 1]] + p.n.car.demand[[i.d]] * p.n.car.return[[i.r]]
      reward[[n.car + 1]] <- reward[[n.car + 1]] + 10.0 * n.car.demand[[i.d]] * p.n.car.demand[[i.d]] * p.n.car.return[[i.r]]
    }
  }
  # print(sum(p.n.car)) # FOR DEBUG: 1でないとおかしい
  return(list(p.n.car=p.n.car, reward=reward))
}

# 方策評価
Evaluate <- function(v, pi, gamma) {
  threshold <- 1e-3 # 適当
  while (TRUE) {
    delta <- 0
    v.org <- v
    for (n1 in 0:20) {
      for (n2 in 0:20) {
        next.s.office1 <- GetNextState(1, n1 - pi[[n1 + 1, n2 + 1]])
        next.s.office2 <- GetNextState(2, n2 + pi[[n1 + 1, n2 + 1]])
        p.dynam <- next.s.office1$p.n.car %o% next.s.office2$p.n.car
        r1.dynam <- next.s.office1$reward %o% next.s.office2$p.n.car
        r2.dynam <- next.s.office1$p.n.car %o% next.s.office2$reward
        r.dynam <- r1.dynam + r2.dynam - 2.0 * abs(pi[[n1 + 1, n2 + 1]]) * p.dynam
        v[[n1 + 1, n2 + 1]] <- sum(r.dynam + p.dynam * gamma * v.org) # Bellman 方程式による更新
        delta <- max(delta, abs(v.org[[n1 + 1, n2 + 1]] - v[[n1 + 1, n2 + 1]]))
      }
    }
    if (delta < threshold) {
      break
    }
  }
  return(v)
}

# 方策改善
Improve <- function(v, pi, gamma) {
  for (n1 in 0:20) {
    for (n2 in 0:20) {
      pi.optimal <- 0
      v.max <- NA
      # この状態でありうる方策を全探索
      pi.min <- - min(5, min(n2, 20 - n1))
      pi.max <- min(5, min(n1, 20 - n2))
      for (pi.cand in pi.min:pi.max) {
        next.s.office1 <- GetNextState(1, n1 - pi.cand)
        next.s.office2 <- GetNextState(2, n2 + pi.cand)
        p.dynam <- next.s.office1$p.n.car %o% next.s.office2$p.n.car
        r1.dynam <- next.s.office1$reward %o% next.s.office2$p.n.car
        r2.dynam <- next.s.office1$p.n.car %o% next.s.office2$reward
        r.dynam <- r1.dynam + r2.dynam - 2.0 * abs(pi.cand) * p.dynam
        v.cand <- sum(r.dynam + p.dynam * gamma * v)
        if (is.na(v.max)) {
          v.max <- v.cand
          pi.optimal <- pi.cand
        } else if (v.cand > v.max) {
          v.max <- v.cand
          pi.optimal <- pi.cand
        }
      }
      # 現在の状態価値関数下でのグリーディ方策に更新
      pi[[n1 + 1, n2 + 1]] <- pi.optimal
    }
  }
  return(pi)
}

# ------------------------------ START MAIN ------------------------------
# 営業終了時に2つの各営業所のレンタカー台数は 0～20台 なので 21x21 個の状態がある
# 状態価値関数と方策 (第1→第2営業所への移送台数、負なら逆向きの移送とする) は0で初期化
v  <- matrix(rep(0, 21*21), nrow=21)
pi <- matrix(rep(0, 21*21), nrow=21)

gamma <- 0.9 # 割引率

while(TRUE) {
  v <- Evaluate(v, pi, gamma) # 方策評価
  pi.org <- pi
  pi <- Improve(v, pi, gamma) # 方策改善
  if (all(pi.org == pi)) { # 方策が更新されなくなったら終わり
    message("FINISHED")
    break
  }
}