反応器計算におけるQ学習

はじめに

簡単な例として、管型反応器において反応温度を調節して所望の反応率を得るという問題について強化学習（Q学習）を適用してみる。

管型反応器の例題

対象とする問題は下記の通り。

• 成分Aを反応して生成物Bを得る管型反応器がある
• 反応器は連続運転しており、1日に1回だけ反応率を確認でき、その際に反応温度を変更することができる
• 反応温度は1回につき0.5℃変更できる（上げることも、下げることも、変えないことも可能）
• 反応温度を変更することで反応率を変化させることができる

ある初期温度で運転を開始したのちに所望の反応率に調整することを学習できるか。

ここで反応率と反応温度の関係は図中の式で表されるものとし、本例題ではk=0.01とする。 反応率と温度は下記グラフの関係となる。

管型反応器のモデル

反応器の状態と行動を次の表の通り整理する。ここでは反応率の目標を0.9599〜0.96としている。 状態は目標より低い、同じ、高い、の3つとし、行動は温度を上げる、何もしない、温度を下げる、の３つとする。 反応率が目標と同じで、かつ何もしないときに報酬が得られるものとする。

import numpy as np
import matplotlib.pyplot as plt

class Tubular(object):
    def __init__(self, temp):
        self._k = 0.01 # kを0.01とする
        ## ログとして[反応温度, 反応率]を残す
        self._log = [[temp, self.calc_conv(self._k, temp)]]

    def calc_conv(self, k, temp):
        ## 反応率を計算
        return 1 - np.exp(-k * temp)

    def state(self):
        ## 状態を返す 反応率が目標内なら1、それより下は0、上は2
        _, conv = self._log[-1] # 一番最新のlogからデータ取得
        if conv < 0.9599:
            st = 0
        elif conv > 0.9600:
            st = 2
        else:
            st = 1
        return st

    def step(self, action):
        ## 計算を1日分進めるステップ計算
        ## まず報酬計算  状態1で行動1(何もしない)のとき報酬あり
        if self.state()==1 and action==1:
            reward = 1
        else:
            reward = 0
        ## 1日分進める
        temp, _  = self._log[-1] # 一番最新のlogからデータ取得
        if action == 0:
            temp += 0.5
        elif action == 2:
            temp -= 0.5
        conv = self.calc_conv(self._k, temp)
        self._log.append([temp, conv])
        return reward

Q値に従った次の行動の決定

def decide_action(qtable, state, episode):
    eps = 0.8 - 0.8*episode/800 # ε-グリーディー法
    if eps <= np.random.uniform(0, 1):
        ## Q値が最大のインデックスを返す
        t = np.where(qtable[state]==qtable[state].max())[0]
    else:
        ## 取りうる行動すべてを返す
        t = np.arange(qtable.shape[-1])
    return np.random.choice(t) # 行動の候補からランダムに選ぶ

Q値の更新

def update_qtable(qtable, state, action, reward, next_state):
    gamma = 0.9
    alpha = 0.5
    next_qmax = max(qtable[next_state])
    qtable[state, action] = (1 - alpha) * qtable[state, action] + alpha * (
        reward + gamma * next_qmax)
    return qtable

学習する

1episodeあたり100日運転することとして、1000episode学習させる。 各episodeの初期の反応温度はランダムに決める。

operation_period = 100 # 100日分
episode_num = 1000
qtable = np.zeros((3, 3)) # サイズは状態数✕行動数

for episode in range(episode_num):
    t = np.random.uniform(300, 350)
    env = Tubular(t)

    for d in range(operation_period):
        ## まず現在の状態を確認
        state = env.state()
        ## 次に行動を決める
        action = decide_action(qtable, state, episode)
        ## 決めた行動に従い運転を1日進める.また報酬を環境から得る
        reward = env.step(action)
        ## 新しい状態を確認
        next_state = env.state()
        ## Qテーブル更新
        qtable = update_qtable(qtable, state, action, reward, next_state)

学習結果

Qテーブルの一例。反応率が目標未満なら温度を上げ、目標より高いなら温度を下げるような値になっている。

#          温度上げる       何もしない      温度下げる
array([[1.97888446e-02, 3.12142611e-12, 1.27102555e-18], # 反応率目標未満
       [2.34887193e+00, 1.00000000e+01, 5.69672696e+00], # 目標内
       [7.11430684e-27, 7.22125479e-27, 1.91366261e-02]]). # 目標より高い

次の結果確認用コードを用いて各種初期温度を計算すると、いずれも最終的に所望の反応率に到達していることを確認できる。

def suisan(init_temp):
    env = Tubular(init_temp)
    for d in range(operation_period):
        state = env.state()
        action = decide_action(qtable, state, 10000) # episode数は適当に大きい値を入れる。探索をしないため。
        _ = env.step(action)
    t, c = np.array(env._log).T
    fig = plt.figure(figsize=(11,4))
    ax1 = fig.add_subplot(121)
    ax2 = fig.add_subplot(122)
    ax1.set_xlabel("Days")
    ax2.set_xlabel("Days")
    ax1.set_ylabel("Temperature")
    ax2.set_ylabel("Conversion")
    ax1.plot(t)
    ax2.plot(c)
    plt.show()

初期温度 310

初期温度 340