ディープラーニングによる化学プラントの異常検知の簡単な例題
はじめに
例えば連続で運転し化学品を製造するプラントにおいて、何らかの原因で反応温度や圧力が異常に上昇したり、製品中の不純物が増加するなどして運転が継続できず、プラント停止を余儀なくされることがある。この場合、プラント設備の復旧や製品在庫管理・販売計画見直しに多大な人手・費用が必要となる。
異常の兆候に早い段階で気付くことができれば、対策を講じることで必要な人手・費用を抑えることができる。このためAI(ディープラーニング)による化学プラントの異常検知が注目されており、例えば下記の三井化学の例がある。ここでは簡単な例題でディープラーニングによる異常検知をやってみる。
人工知能（AI）を用いて、化学プラントの製造過程で製品の品質予測に成功｜2016｜ニュースリリース｜三井化学株式会社
例題
下記の原料Aから生成物Bを製造する反応器を考える。①原料流量、②原料温度、③反応圧力、④熱媒体温度を制御して、⑤出口温度をある一定の温度範囲内に収めるように運転している（範囲から外れると製品品質悪化)。
以下のグラフはある1000日間の運転データである。900日を過ぎたあたりで出口温度がこのプラントの管理上限値である320℃を超え異常が認識された。出口温度の上昇は反応副生物の分解によるもので、急激に反応が起こるのでなく、ゆっくりと時間をかけて分解量が増加していくことがわかっているものとする。
反応量が少ないうちに出口温度の僅かな上昇で異常を検知できることが好ましいが、制御値①②③④の振れが出口温度を振れさせてしまい、副生物分解による僅かな温度上昇を見分けることは難しい状況である。
ディープラーニングを用いて早い段階で異常を検知することができるか？
※このデータは下記コードから人工的に生成した。気になる人は参考に。制御値①②③④はランダムで振れを発生させ、出口温度は①②③④を変数とした所定の関数から算出している。さらに700日目から日数に応じて増加する異常値を出口温度に加えている。ここで使用したデータそのものもリンク先に置いてある。
異常検知の手順
次の手順でディープラーニングを用いて異常を検知する。
- 正常な期間のデータをディープラーニングで学習し、制御値①②③④と出口温度を相関するモデルを作る
- 上記1.のモデルを使って異常な期間の出口温度の予測値を得る
- 異常な期間の実際の出口温度とディープラーニングで予測した出口温度の差を異常度とする
正常期間データを学習したモデルから得た異常期間の予測値は副生物分解による温度上昇を含まない値となる。したがって、この値と異常期間の実際の値の差を取ることで異常さの度合いが見えるはずである。
正常期間データの学習
700日目までは正常だったと考えて、その期間のデータを再現するようネットワークを学習させた。
なお用いたネットワークは下記の通りで、ノード数20の中間層を2つ入れた。
## ニューラルネットワークモデル設定 class Model(nn.Module): def __init__(self): super(Model, self).__init__() self.l1 = nn.Linear(4, 20) # 入力4つ (原料流量、原料温度、反応圧力、熱媒体温度) self.l2 = nn.Linear(20, 20) self.l3 = nn.Linear(20, 1) # 出力1つ (出口温度) def forward(self, x): x = F.relu(self.l1(x)) x = F.relu(self.l2(x)) return self.l3(x)
学習させた結果、次のグラフとなった。色々と試行錯誤したものの初期の乖離は解消されなかった。しかしながら、それ以外の期間についてはほぼデータを再現できている。
学習したモデルによる出口温度予測および異常度の算出
異常と考えられる期間も含めたデータとモデル予測値は次のグラフとなる。異常を認識したタイミングよりも早い時点でデータとモデル予測値に差が生じていることが見て取れる。
データとモデル予測値の差の絶対値を「異常度」と定義し、それをグラフにプロットすると次の図となる。異常度がいくつになったら異常と判断するかがポイントとなってくるが、例えば初期を除き10℃以上を異常と考えるなら、100日程度早めに異常兆候を知ることができる。100日あれば事前に十分な対策を講じることができるだろう。
おわりに
実際のプラントでは様々なノイズが制御値や計測値に影響するため、こう簡単にはいかないかもしれないが、実務でも適用をトライしたいと思う。
使用したコード
深層強化学習の簡単な例 〜Double DQN適用〜
はじめに
以前の記事で簡単な問題を対象にDQN、Experience Replayを適用してみたが、今回はさらにDouble DQNを追加してみる。なお下記サイトを参考にした。
Let’s make a DQN: Double Learning and Prioritized Experience Replay – ヤロミル
https://www.renom.jp/ja/notebooks/product/renom_rl/ddqn/notebook.html
Double DQN
通常のDQNでは次の形でニューラルネットワークを更新している。
ポイントはで、次の状態の最大のQ値を求めるのに自分自身のネットワークを使用している点である。もし仮にQ値の初期値が変に偏っていて、ある行動のQ値が他に比べ極端に高い場合、同じネットワークで更新を繰り返すと、その行動の価値が過大に評価されてしまうことになる。
これを緩和するために提案された方法がDouble DQNである。この方法では、メインのネットワークとは別のネットワーク（Target network)を作り、ネットワーク更新時に次の状態における行動はメインネットワークで決めるが、その行動の価値評価は別のネットワークで行う。なお別のネットワークには過去のメインネットワークを用いる。この方法を用いる場合のネットワーク更新式は次の形となる。
対象とする問題
前回記事と同一の報酬払出装置を考える。前回のプログラムを修正してDouble DQNを適用する。
深層強化学習(Deep Q Network, DQN)の簡単な例 〜Experience Replay追加〜 - 化学系エンジニアがAIを学ぶ
Double DQN の適用
次の通り、メインネットワークに加え、メインとは別のネットワーク(Target network)を準備する。
import torch import torch.nn as nn import torch.nn.functional as F class DQN(nn.Module): def __init__(self): super(DQN, self).__init__() self.l1 = nn.Linear(1, 3) self.l2 = nn.Linear(3, 3) self.l3 = nn.Linear(3, 2) def forward(self, x): x = F.relu(self.l1(x)) x = F.relu(self.l2(x)) x = self.l3(x) return x dqn = DQN() # メインネットワーク dqn_t = DQN() # もう１つのネットワーク(Target network) optimizer = torch.optim.SGD(dqn.parameters(), lr=0.01) criterion = nn.MSELoss()
Double DQNを適用する場合のネットワーク更新は下記の形となる。通常のDQNと異なるのは下記「Q値の算出」にあるa_m、max_q_nextを求める部分だけである。（なおここではExperience Replayを利用している。）
BATCH_SIZE = 50 def update_dqn(replay_memory): ## メモリがバッチサイズより小さいときは何もしない if len(replay_memory) < BATCH_SIZE: return ## ミニバッチ取得 transitions = replay_memory.sample(BATCH_SIZE) ## (状態、行動、次の状態、報酬)✕バッチサイズ を (状態xバッチサイズ、行動✕バッチサイズ、、、)に変える batch = Transition(*zip(*transitions)) ## 各値をtensorに変換 state = torch.FloatTensor(batch.state).unsqueeze(1) action = torch.LongTensor(batch.action).unsqueeze(1) next_state = torch.FloatTensor(batch.next_state).unsqueeze(1) reward = torch.FloatTensor(batch.reward).unsqueeze(1) ## Q値の算出 q_now = dqn(state).gather(-1, action) # 今の状態のQ値 a_m = dqn(next_state).max(-1)[1].unsqueeze(1) # メインネットワークで決めた次のaction,[1]でindexを得る max_q_next = dqn_t(next_state).gather(-1, a_m) # Target networkから上で決めたactionのQ値を得る gamma = 0.9 q_target = reward + gamma * max_q_next # 目標のQ値 # DQNパラメータ更新 optimizer.zero_grad() loss = criterion(q_now, q_target) # 今のQ値と目標のQ値で誤差を取る loss.backward() optimizer.step() return loss.item() # lossの確認のために返す
学習時に所定の計算ステップごとにTarget networkにメインネットワークをコピーして更新される形とする（下記コードの一番下）。更新されるまでは学習時に過去のメインネットワークを参照する形となる。
NUM_EPISODES = 1200 NUM_STEPS = 5 log = [] # 結果のプロット用 for episode in range(NUM_EPISODES): env = Dispenser(0) total_reward = 0 # 1エピソードでの報酬の合計を保持する for s in range(NUM_STEPS): ## 現在の状態を確認 state = env.state ## 行動を決める action = decide_action(state, episode) ## 決めた行動に従いステップを進める。次の状態、報酬を得る next_state, reward = env.step(action) ## 結果をReplayMemoryに記録 memory.push(state, action, next_state, reward) ## DQNを更新 loss = update_dqn(memory) total_reward += reward log.append([total_reward, loss]) if (episode+1)%100==0: print(episode+1, loss) ## 2エピソード = 計算10ステップ 毎の頻度でTarget networkにメインネットワークをコピー if episode%2==0: dqn_t.load_state_dict(dqn.state_dict())
学習結果
多少計算が安定している感じがするが、このくらいの問題ではあまり差はないのかもしれない。
コード
メモ: PyTorch 重み、バイアスの初期化について
はじめに
PyTorchのニューラルネットワークの重み・バイアスの初期化についてのメモを記す。
重み
重みの内容は次のようにして確認できる。
>>> import torch.nn as nn >>> l = nn.Linear(1, 3) >>> l.weight Parameter containing: tensor([[ 0.6204], [-0.5651], [-0.6809]], requires_grad=True)
重みの初期化は次のようにnn.initモジュールの関数を使用することで実施できる。
>>> nn.init.xavier_uniform_(l.weight) # Xavierの初期値 Parameter containing: tensor([[ 0.5318], [-0.8438], [-1.1093]], requires_grad=True) >>> nn.init.kaiming_uniform_(l.weight) # Heの初期値 Parameter containing: tensor([[ 2.3349], [-0.9218], [-1.3357]], requires_grad=True)
単純な正規分布・一様分布や定数を設定できる関数もnn.init内に準備されている。
詳しくはこちらを参照: https://pytorch.org/docs/master/nn.init.html
バイアス
バイアスの内容は重みと同様に確認できる。
>>> l.bias Parameter containing: tensor([-0.9920, 0.6406, -0.6558], requires_grad=True)
バイアスの初期化も重みと同様のやり方。
>>> nn.init.uniform_(l.bias, 0, 1) # 一様分布 Parameter containing: tensor([0.8032, 0.7810, 0.1619], requires_grad=True) >>> nn.init.constant_(l.bias, 0) # 定数 Parameter containing: tensor([0., 0., 0.], requires_grad=True)
デフォルトの初期値
nn.Linearのソースを見ると重みについてはHeの初期値が入っている模様。
pytorch/linear.py at master · pytorch/pytorch · GitHub
def reset_parameters(self): init.kaiming_uniform_(self.weight, a=math.sqrt(5)) if self.bias is not None: fan_in, _ = init._calculate_fan_in_and_fan_out(self.weight) bound = 1 / math.sqrt(fan_in) init.uniform_(self.bias, -bound, bound)
参考
python - How to initialize weights in PyTorch? - Stack Overflow
What is the default initialization of a conv2d layer and linear layer? - PyTorch Forums
ニュース: Preferred Networks、JXTGエネルギーと石油精製プラントの最適化・自動化に関する共同研究を開始
Preferred NetworksがJXTGエネルギーと石油精製プラントの最適化・自動化に関する共同研究を開始したそうです。１０億円もの資金が調達されるそうで、一気に技術開発が進むことが期待されます。近い業界に関わる者としては、得られた知見などを可能な範囲でオープンにしてもらえたらと思いますが、難しいでしょうね。。。
深層強化学習(Deep Q Network, DQN)の簡単な例 〜Experience Replay追加〜
はじめに
前回の記事でOpenAI Gymを使わず非常に簡単な問題を対象にDQNを適用してみたが、"Experience Replay"を入れていなかった。今回は前回の問題にExperience Replayを追加してみる。なおこれを実施するにあたり、下記サイトを参考にした。
第15回 CartPole課題で深層強化学習DQNを実装｜Tech Book Zone Manatee
Reinforcement Learning (DQN) Tutorial — PyTorch Tutorials 1.5.1 documentation
Experience Replay
行動とその結果を経験として記録しておき、その過去の経験をサンプリングして学習データとする方法。通常の逐次で学習する方法に比べ、時間的に離れたデータを使用することができ、偏りを抑えた安定した学習ができるとのこと。複数データをまとめてニューラルネットワークにミニバッチ学習させることもできる。
具体的には「現在の状態」、「行動」、「行動後の状態」、「報酬」のデータの組をある一定の数だけメモリに残しておいて、その中から学習に使うデータをいくつかランダムにサンプリングし、これらを１つのバッチとしてニューラルネットワークに学習させる。
対象とする問題
前回記事と同一の報酬払出装置を考える。プログラムは前回のものにいくつか追加・修正してExperience Replayを適用する。
深層強化学習(Deep Q Network, DQN)の簡単な例 - 化学系エンジニアがAIを学ぶ
Experience Replayの実装
下記のpytorch.orgのチュートリアルにあるように、namedtupleを用いて経験を記録するクラスを追加する。
Reinforcement Learning (DQN) Tutorial — PyTorch Tutorials 1.5.1 documentation
ここでは500回分の行動の結果を記録することにしている(数は適当)。500を超えると古い方から順番に上書きされる形となっている。サンプリング時は指定したバッチサイズ数だけのデータが取り出される。
from collections import namedtuple # 「現在の状態」、「行動」、「行動後の状態」、「報酬」を記録するnamedtuple Transition = namedtuple('Transition', ('state', 'action', 'next_state', 'reward')) import random class ReplayMemory(object): def __init__(self, capacity): self.capacity = capacity # メモリの許容サイズ self.memory = [] # 経験を保存するリスト self.index = 0 # 保存した経験リストのindex def push(self, state, action, next_state, reward): # 経験をリストに追加(push)する if len(self.memory) < self.capacity: # 許容サイズ以下なら追加 self.memory.append(None) self.memory[self.index] = Transition(state, action, next_state, reward) self.index = (self.index + 1) % self.capacity # indexを1つ進める。許容サイズ超えたら古い順から上書き def sample(self, batch_size): # バッチサイズ分の経験のランダムサンプリング return random.sample(self.memory, batch_size) def __len__(self): return len(self.memory) CAPACITY = 500 memory = ReplayMemory(CAPACITY)
DQNの学習
上記の行動を記録したクラスを利用できる形に、DQNの学習部分を修正する。
BATCH_SIZE = 50 # サンプリングするデータ数 def update_dqn(replay_memory): # 行動を記録したクラスを引数とする ## メモリがバッチサイズより小さいときはまだ学習しない if len(replay_memory) < BATCH_SIZE: return ## バッチ取得 transitions = replay_memory.sample(BATCH_SIZE) ## (状態、行動、次の状態、報酬)✕バッチサイズ を (状態xバッチサイズ、行動✕バッチサイズ、、、)に変える batch = Transition(*zip(*transitions)) ## 各値をtensorに変換 state = torch.FloatTensor(batch.state).unsqueeze(1) action = torch.LongTensor(batch.action).unsqueeze(1) next_state = torch.FloatTensor(batch.next_state).unsqueeze(1) reward = torch.FloatTensor(batch.reward).unsqueeze(1) ## Q値の算出 q_now = dqn(state).gather(-1, action) # 今の状態のQ値 max_q_next = dqn(next_state).max(-1)[0].unsqueeze(1).detach() # 状態移行後の最大のQ値 gamma = 0.9 q_target = reward + gamma * max_q_next # 目標のQ値 # DQNパラメータ更新 optimizer.zero_grad() loss = criterion(q_now, q_target) # 今のQ値と目標のQ値で誤差を取る loss.backward() optimizer.step() return loss.item() # lossの確認のために返す
学習結果
他のコードは前回記事と同じである。下図に学習の推移を示す。左は誤差の推移、右は1エピソードにおけるトータル報酬の推移である。実際に計算してみるとわかるが逐次で学習させるより、Experience Replayを適用するほうが学習が安定しやすい。
コード
参考
・第15回 CartPole課題で深層強化学習DQNを実装｜Tech Book Zone Manatee
・Reinforcement Learning (DQN) Tutorial — PyTorch Tutorials 1.5.1 documentation
深層強化学習(Deep Q Network, DQN)の簡単な例
はじめに
DQNを学ぼうとして色々と調べたが、どこもかしこもOpenAI Gymを使っていて、まずそれの扱いから考えないといけないのでつらい。ここではOpenAI Gymを使わず、非常に簡単な問題を対象にDQNを適用してみることとする。Python、PyTorchを用いる。
また、DQNをやろうとすると出てくる"Experience Replay"。これも入れようとすると考えるのがつらいので、入れない。
下記記事にDQN全般やExperience Replayについてちょっと解説がある。
Deep Q Network (DQN) - DeepLearningを勉強する人
対象とする問題
下図のような報酬払出装置を考える。装置には「電源」ボタンと「払出」ボタンが設置されている。「電源」ボタンを押すと装置電源のON/OFFが切り替わる。電源ONのときに「払出」ボタンを押すと報酬が払い出される。電源OFFのときに「払出」ボタンを押しても報酬は払い出されない。
1回の行動で、いずれかのボタンを1度押すことができるものとする。5回の行動を行う場合に、報酬を最大化する手順を学習することができるか？
なお、この問題の内容に関して下記書籍を参考にした。
この払出装置(Dispenser)をプログラムにすると次の通りとなる。
class Dispenser(object): def __init__(self, init_state): """ 初期のON/OFF状態を設定する init_state: 0->電源OFF、1->電源ON """ self.state = init_state def powerbutton(self): """ 電源ボタンを押し、ON/OFFを切り替える """ if self.state == 0: self.state = 1 else: self.state = 0 def step(self, action): """ 払出機を操作する action: 0->電源ボタンを押す 1->払出ボタンを押す 状態と報酬が返る """ if action == 0: # 電源ボタンを押した場合 self.powerbutton() # 電源ON/OFF切り替え reward = 0 # 報酬はない else: # 払出ボタンを押した場合 if self.state == 1: reward = 1 # 電源ONのときのみ報酬あり else: reward = 0 return self.state, reward
Deep Q Network (DQN)
報酬払出機の状態は2つ、取れる行動は2つであるので、これを表すQテーブルは次のようになる。
Q学習ではこのQテーブルのQ値（Q(0,0), Q(0,1), Q(1,0), Q(1,1)）を逐次更新していくが、DQNではこれらのQ値を得る関数を逐次更新していく。その関数を表すのにニューラルネットワークを用いるが、次のようなイメージとなる。
状態を入力として、Q値を出力として得る形である。ここでは上図の通りのニューラルネットワークを用いることとする。（あまりディープではないが・・・）
import torch import torch.nn as nn import torch.nn.functional as F class DQN(nn.Module): def __init__(self): super(DQN, self).__init__() self.l1 = nn.Linear(1, 3) self.l2 = nn.Linear(3, 3) self.l3 = nn.Linear(3, 2) def forward(self, x): x = F.relu(self.l1(x)) x = F.relu(self.l2(x)) x = self.l3(x) return x dqn = DQN() optimizer = torch.optim.SGD(dqn.parameters(), lr=0.01) criterion = nn.MSELoss()
DQNの学習方法
では、どのような情報を使ってこのニューラルネットワークを学習させていくのか、ということになるが、これにはQ学習のQ値の更新式を用いる。
この式は今のQ値（）を、今の推定の目標値（
）に近づくよう更新するものである。すなわち、
となることを目指している。よってと
の誤差を小さくするようにニューラルネットワークを学習させるという形にすればよいことになる。いずれの値も現在の状態、行動、次の状態、報酬が分かれば算出可能であり、DQNの学習をプログラムにすると次のようになる。
def update_dqn(state, action, next_state, reward): ## 各変数をtensorに変換 state = torch.FloatTensor([state]) action = torch.LongTensor([action]) # indexとして使うのでLong型 next_state = torch.FloatTensor([next_state]) ## Q値の算出 q_now = dqn(state).gather(-1, action) # 今の状態のQ値 max_q_next = dqn(next_state).max(-1)[0].detach() # 状態移行後の最大のQ値 gamma = 0.9 q_target = reward + gamma * max_q_next # 目標のQ値 # DQNパラメータ更新 optimizer.zero_grad() loss = criterion(q_now, q_target) # 今のQ値と目標のQ値で誤差を取る loss.backward() optimizer.step() return loss.item() # lossの確認のために返す
DQNからの行動の決定
DQNは入力された状態に対し、取りうるすべての行動のQ値のtensorを返す（1次元tensor)。そこから最大のQ値のindexを取り出して(0または1)、それをそのまま行動としている。
import numpy as np EPS_START = 0.9 EPS_END = 0.0 EPS_DECAY = 200 def decide_action(state, episode): state = torch.FloatTensor([state]) # 状態を1次元tensorに変換 ## ε-グリーディー法 eps = EPS_END + (EPS_START - EPS_END) * np.exp(-episode / EPS_DECAY) if eps <= np.random.uniform(0, 1): with torch.no_grad(): action = dqn(state).max(-1)[1] # Q値が最大のindexが得られる action = action.item() # 0次元tensorを通常の数値に変換 else: num_actions = len(dqn(state)) # action数取得 action = np.random.choice(np.arange(num_actions)) # ランダム行動 return action
学習する
現在の状態からDQNより行動を決めて、環境から次の状態・報酬を取得し、これらの結果を用いてDQNを更新する、を繰り返す。エピソードの数や上記のoptimizerの学習率、ε-グリーディー法のepsの変化のさせ方はうまく学習が進むように調節して決めている。（適度な値を決めるのにかなり苦労した。)
import matplotlib.pyplot as plt NUM_EPISODES = 1200 NUM_STEPS = 5 log = [] # 結果のプロット用 for episode in range(NUM_EPISODES): env = Dispenser(0) total_reward = 0 # 1エピソードでの報酬の合計を保持する for s in range(NUM_STEPS): ## 現在の状態を確認 state = env.state ## 行動を決める action = decide_action(state, episode) ## 決めた行動に従いステップを進める。次の状態、報酬を得る next_state, reward = env.step(action) ## DQNを更新 loss = update_dqn(state, action, next_state, reward) total_reward += reward log.append([total_reward, loss]) ## 結果表示 r, l = np.array(log).T fig = plt.figure(figsize=(11,4)) ax1 = fig.add_subplot(121) ax2 = fig.add_subplot(122) ax1.set_xlabel("Episode") ax1.set_ylabel("Loss") ax1.set_yscale("log") ax2.set_xlabel("Episode") ax2.set_ylabel("Total reward") ax1.plot(l) ax2.plot(r) plt.show()
学習結果
下図に学習の推移を示す。左は誤差の推移、右は1エピソードにおけるトータル報酬の推移である。払出機は初め電源OFFのため、5回の行動で得られる最大のトータル報酬は4であるが、右下のグラフから最終的にトータル報酬4が得られる行動を適切に学習できていることを確認できる。
学習後のDQNから各状態におけるQ値を確認できるが、今回の学習での値は下記である。
>>> dqn(torch.FloatTensor([0])) tensor([9.0000, 8.1200], grad_fn=<AddBackward0>) >>> dqn(torch.FloatTensor([1])) tensor([8.1060, 10.0000], grad_fn=<AddBackward0>)
計算の詳細は示さないが、今回と同一の問題にQ学習を適用するとこれらとほぼ同じQ値が得られ、DQNによってQ値を適切に表現する関数を構築できていることを確認できる。
さいごに
DQNを適用するのが全く意味がないような簡単な問題にDQNを適用した。"Experience Replay"も考えなかった（状態✕行動の組み合わせが4パターンしかない問題なのであまり効果はないのかも）。簡単すぎる問題であるがDQNを適用する上での基礎的なポイントは掴めたのではないかと思う。
コード
ここで示した全コードおよびQ学習を適用した場合のコードはこちら:
参考
・Reinforcement Learning (DQN) Tutorial — PyTorch Tutorials 1.1.0 documentation
反応器計算におけるQ学習2
はじめに
以前の記事で管型反応器において所望の反応率を得るという問題について強化学習（Q学習）を適用したが、 今回は反応率に応じて目的成分の選択率が変化するという条件で、収率を最大化する問題を対象とする。
管型反応器の例題
対象とする問題は下記の通り。
- 成分Aを反応して目的の生成物Bを得る管型反応器がある
- 反応器は連続運転しており、1日に1回だけAの反応率とBの収率を確認でき、その際に反応温度を変更することができる
- 反応温度は1回につき0.5℃変更できる（上げることも、下げることも、変えないことも可能）
- 反応温度を変更することで反応率を変化させることができる
- 生成物Bと一緒に副生物Cも生成し、反応率に応じて生成物Bの選択率は変化する
ある初期温度で運転を開始したのちに、収率を最大化することを学習できるか?
ここで反応温度・反応率・選択率・収率の関係は図中の式で表されるものとし、 本例題ではk=0.01、a=100とする。 反応温度と反応率、反応率と選択率の関係は下記グラフとなる。
なお、収率が最大となる条件はこの時点で計算可能であり、反応率0.955、収率0.945である。 これを強化学習で推算できるかがポイントとなる。
管型反応器のモデル
反応器の状態と行動を次の表の通り整理する。 状態は反応率を離散化して表現する。ここでは反応率0.9未満を状態0とし、0.9から1.0を100分割して番号を振った。 (収率が最大となる反応率が少なくとも0.9以上とわかっている前提）
報酬はとし、1に満たない分をマイナスの報酬（ペナルティ）として与える形とした。
import numpy as np import matplotlib.pyplot as plt class Tubular_selec(object): def __init__(self, temp): self._k = 0.01 # kを0.01とする self._a = 100.0 # aを100とする init_conv = self.calc_conv(self._k, temp) init_yield = self.calc_yield(self._a, init_conv) ## ログとして[温度, 反応率, 収率]を残す self._log = [[temp, init_conv, init_yield]] def calc_conv(self, k, temp): ## 反応率計算 return 1 - np.exp(-k * temp) def calc_selc(self, a, conv): ## 選択率計算 return 1 - conv**a def calc_yield(self, a, conv): ## 収率計算 return conv * self.calc_selc(a, conv) def state(self): ## 状態を返す 状態を転化率0.9-1.0の間で100分割する ## 0.9未満を状態0とするので、状態は全部で101個 _, conv, _ = self._log[-1] # 一番最新のlogからデータ取得 n = 100 if conv < 0.9: # 反応率 0.1未満は状態0とする st = 0 else: st = int(np.ceil((conv - 0.9)*10*n)) # 0.9~1.0をn分割 return st def step(self, action): temp, _, _ = self._log[-1] # 最新の温度を取得 ## 行動0で温度UP、行動2で温度Down、行動1は何もしない if action == 0: temp += 0.5 elif action == 2: temp -= 0.5 conv = self.calc_conv(self._k, temp) yild = self.calc_yield(self._a, conv) self._log.append([temp, conv, yild]) # 1.0に満たない分を負の報酬(ペナルティ)として与える reward = - abs(1.0 - yild) return reward
Q値に従った次の行動の決定
def decide_action(qtable, state, episode): eps = 1 - episode/240 # ε-グリーディー法 if eps <= np.random.uniform(0, 1): # グリーディー行動 t = np.where(qtable[state]==qtable[state].max())[0] # Q値が最大のインデックスを返す 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日運転することとして、300episode学習させる。
OPERATION_PERIOD = 100 # 100日分 EPISODE_NUM = 300 qtable = np.zeros((101, 3)) # サイズは状態数✕行動数 ## 学習状況確認用グラフ設定 fig = plt.figure(figsize=(11,4)) ax1 = fig.add_subplot(121) ax2 = fig.add_subplot(122) ## 学習 for episode in range(EPISODE_NUM): #t = np.random.uniform(280, 350) # 初期温度をランダムにする場合 env = Tubular_selec(280) 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) ## 10episodeにつき1回のグラフ描画 if (episode+1)%10 == 0: t, c, y = np.array(env._log).T ax1.cla() ax2.cla() ax1.set_xlabel("Days") ax1.set_ylabel("Temperature") ax1.set_ylim(270,320) ax2.set_xlabel("Days") ax2.set_ylabel("Conversion, Yield") ax2.set_ylim(0.93,0.96) ax1.plot(t) ax2.plot(c, label="Conversion") ax2.plot(y, label="Yield") ax2.legend() plt.pause(0.1)
学習結果
収率最大となるのは反応率0.955、収率0.945のときであるが、 反応温度を上げて、まっすぐその状態に近づいていることが確認できる。
