連続晶析における粒度分布遷移過程の状態空間モデルによる分析

(1)

＜修士論文＞

連続晶析における粒度分布遷移過程の

状態空間モデルによる分析

滋賀大学大学院

データサイエンス研究科

データサイエンス専攻

修了年度：２０２０年度

学籍番号：６０１９１０１

氏名：秋山浩希

指導教員：河本薫

提出年月日：２０２１年

1 月２０日

(2)

第1章はじめに 3 1.1 背景と目的 . . . 3 1.2 本論文の構成 . . . 4 第2章晶析と実験装置 5 2.1 晶析 . . . 5 2.2 実験装置 . . . 9 2.3 運転プログラム . . . 14 2.4 粒径の計測方法 . . . 15 第3章実験条件と実験結果 17 3.1 実験条件 . . . 17 3.2 実験結果 . . . 18 第4章状態空間モデルによる分析 21 4.1 探索的データ解析 . . . 21 4.2 状態空間モデルの構築 . . . 26 4.3 ベイズ法によるパラメータ推定 . . . 30

(3)

第5章分析結果と考察 34 5.1 MCMC法による計算結果の収束診断 . . . 34 5.2 推定されたパラメータによるモデルと実験データの比較 . . . 36 5.3 考察 . . . 39 第6章今後の課題と結論 40 6.1 今後の課題 . . . 40 6.2 結論 . . . 43 謝辞 44 参考文献 45

(4)

第

₁

章

はじめに

1.1 背景と目的

本研究は、滋賀大学と住友金属鉱山株式会社（以下「住友金属鉱山」）の共同研究である。住友金属鉱山では、晶析と呼ばれる技術を用いて様々な結晶を製造し、顧客へと納品している。納品する結晶は、その品質を担保しなければならない。すなわち、要求される品質の結晶を製造できなければそれは不良品となり、損失を被ることになる。製造された結晶の品質は、結晶の粒度などによって管理されている。結晶の粒度とは、結晶の大きさやそのばらつき具合のことで、粒度分布を用いて表現される。製造現場では、この粒度分布について、品質面で許容できるかどうかを判断するために、粒度分布の統計量を管理指標としている。統計量と言うと、平均値や標準偏差などが代表的であるが、晶析においては、D50という独特の統計量を用いることが多い。D50の定義については、2 章にて詳しく述べる。一方、晶析をおこなうための晶析装置の運転条件が粒度に及ぼす影響については解明されていないことが多い。実際、製造現場においては、D50を目標レンジに近づけるために運転条件をどのように変えればいいかは暗中模索的であり、目標レンジを超過してしまったり、目標レンジに近づきつつあることに気付かずに運転条件を更に変更してしまうなど、製造準備に長時間を要してしまい、生産性を下げている。このような現場課題に対して、本研究では、D50という統計量ではなく粒度分布を分析対象とする方針をとる。なぜなら、粒度分布を管理する方法がわかればD50の管理もできるようになると考えられるためである。したがって本研究の目的は、晶析装置の運転条

(5)

件が粒度分布にどのような影響を与えるかを分析することによって、結晶品質の安定化方法、すなわち粒度の安定化方法を探ることとした。なお、本研究に用いる晶析装置は、住友金属鉱山より提供いただいた晶析の実験装置を使用する。

1.2 本論文の構成

第2章では、はじめに晶析とは何かについて説明し、つぎに本研究において使用する晶析の実験装置について説明する。第3章では、本研究で用いた実験データを取得する際の実験条件と、その実験結果について説明する。第4章では、実験により得られたデータから探索的データ解析による仮説の構築をおこない、その仮説を状態空間モデルによって表現する。モデルのパラメータはベイズ法を用い、MCMC法によって推定する。第5章では、MCMC法の計算結果の収束診断をおこない、つづいて推定されたパラメータを用いたモデルによる再現データと実験データの比較をおこなう。第6章では、今後の課題として、本研究において未検証となってしまったことや、新たに発見された課題などについて述べる。

(6)

第

₂

章

晶析と実験装置

本章では、晶析とはなにかについて、および本研究において使用する実験装置について説明する。なお本章で説明する内容は、書籍「晶析工学」[1]やオンライン上の解説記事 [2] [3] [4]、および住友金属鉱山の佐藤健司氏よりご教授いただいた知識などにもとづいている。

2.1 晶析

晶析とは、我々が望む品質や特性を持つ結晶を、経済的かつ大量に生産するための技術である。結晶とは、「化学物質の最小単位である原子や分子が、規則正しく積み重なってできた固体」[4]である。たとえば、食塩や砂糖などの一粒一粒が結晶である。このような結晶の特性は、個々の結晶の大きさや、集合体の大きさのばらつきなどの性質によって決まる。たとえば砂糖の結晶であれば、中双糖（ザラメ）のように粗い粒は口内で溶けにくくコクのある味わいとなり、上白糖のように細かい粒は口内で溶けやすく強い甘みを感じる味わいとなる。中双糖と上白糖を混合したような砂糖は別の味わいになるだろうし、混合の割合が異なれば味も異なる。すなわち、同一物質の結晶でも、結晶の大きさやそのばらつきなどが異なる場合、異なる特性を持ちえる。このような理由から、結晶の大きさやそのばらつきを結晶の粒度（粒の度合い）とすれば、結晶の粒度を管理することが、結晶の品質管理では重要になるのである。本節では、まず、結晶の粒度を見るための粒度分布について説明する。つぎに、製造現場で用いられる晶析方法ついて説明する。

(7)

2.1.1 結晶の粒度分布

粒度分布とは、「測定対象となるサンプル粒子群の中に、『どのような大きさ（粒子径）の粒子が、どのような割合（全体を100％とする相対粒子量）で含まれているか』を示す指標（表現手段）」[3]である。粒度分布は、一義的なものではない。まず、結晶は立体的な粒であるため、その大きさ（すなわち粒度）を計測する尺度は、長さや表面積、体積など多義的である。また、割合を計算する場合も、該当する結晶群の個数が全結晶の個数に占める割合（以下、個数割合と呼ぶ）を計算するやり方だけでなく、該当する結晶群の体積が全結晶の体積に占める割合（以下、体積割合と呼ぶ）を計算するやり方もある。その中でも、製造現場における品質管理でよく用いられる粒度分布を図2.1に示す。図2.1: 粒度分布の例以下では図2.1の横軸と縦軸について説明する。図2.1右図のD50については、次項で説明する。はじめに、図2.1の横軸である粒径について説明する。粒径とは、結晶の形状が（真）球であると仮定したときの、球の直径である。ほとんどの粒子の形状は、球や立方体といった単純かつ定量的に表現できるものではなく、複雑かつ不規則であり、直接的に粒子径を定義することはできない。そこで、「球相当径」という間接的な定義を用いる[3]。なおこ

(8)

こでは横軸に球相当径の直径である粒径を用いているが、結晶がもつ特性によっては球相当径の表面積や体積を用いる場合もある。つぎに図2.1の縦軸について説明する。粒径の頻度分布（左図）の縦軸は、すべての結晶のうち、一定の粒径区分（ビンの幅）に属する結晶の割合を示している。体積の累積分布（右図）の縦軸は、「すべての結晶の体積の総和」に対する「ある粒径以下の結晶の体積の総和」の累積相対割合である。この累積分布においては球相当径の体積が用いられることが多い。なお先ほどと同様に、結晶がもつ特性によっては個数や表面積の累積分布で表現する場合もある。

2.1.2 粒度管理指標「

D50

」の定義

D50とは、体積累積分布の中央値の粒径である（図2.1右図）。すなわち、「粒径がD50 以下の結晶の体積の総和」が「すべての結晶の体積の総和の50%」になるということである。D10やD90といった統計値についても、D50と同様に、その粒径以下の結晶の体積の総和がすべての結晶の体積の総和の10%や90%となる粒径として定義される。本論文では、粒度の管理指標として、D50だけを取り上げる。

2.1.3 晶析の方法

何の結晶を生成するか、あるいはどのような特性をもつ結晶を生成するかなどによって、晶析の方法は異なる。ここでは晶析の方法の違いについて、（１）過飽和生成法の違いと（２）装置形式の違い、これら2つのレベルで分類して説明する。（１）過飽和生成法過飽和生成法の違いとは、晶析媒体となる溶液（融液や気相の場合もある）を、どのような方法を用いて過飽和状態にするかという違いである。過飽和状態とは何かについて説明するために、まずは溶解度および飽和状態とは何かについて説明する。溶解度とは、一定の条件下において溶媒100gあたりに溶ける溶質の量である。飽和状態とは、一定の条件下において限界まで溶質を溶液に溶かし、溶質と溶液が平衡状態になっている状態である。また、飽和状態の溶液の濃度を、飽和濃度という。過飽和状態とは、溶解度以上に溶質を溶液に溶かした不安定な状態であり、何らかの衝撃が加わると結晶の析出などが起こる状態である。すなわち、晶析によって結晶を析出させるためには溶液を過飽和状態にしな

(9)

表2.1: 過飽和生成法の表名称方法冷却法溶解度が温度とともに増加する場合には、温度を下げると飽和濃度が低下する。したがって、冷却によって飽和状態とすることができる蒸発法蒸発によって溶液を濃縮することによって過飽和状態を作る貧溶媒添加法良溶媒溶液に貧溶媒を添加して溶解度を下げ、過飽和状態を作る。貧溶媒*に良溶媒溶液を添加する場合もある pH-シフト法アミノ酸水溶液のようにpHによって溶解度が変化する場合には、pHを変化させることにより溶解度を下げ、過飽和状態とすることができる化学反応法化学反応により目的物質を生成させることにより、過飽和状態を作る *溶質を溶かしにくい溶媒を貧溶媒という　　　　　　（「晶析工学」, p.25）ければならず、過飽和状態を作りだすための反応方法の違いを過飽和生成法の違いと呼んでいる。過飽和状態を作りだすための反応方法には様々な方法があるが、ここでは表2.1に主な5つの方法とその名称をまとめた。なお本研究においては、化学反応法という過飽和生成法を用いて晶析をおこなっている。化学反応法とは、化学反応により目的物質を生成させることで過飽和状態を作る過飽和生成法である。本研究において扱う化学物質に関しては、次節にて後述する。（２）装置形式装置形式の違いとは、晶析をおこなう晶析装置（操作）形式の違いのことであり、主に連続式と回分式の2つの形式がある[5]。連続式と回分式の両面を取り入れた半回分式という形式もあるが、特殊な場合や実用されない場合が多いためここでは触れない[6]。連続式は、溶液を連続的に晶析装置に供給し、製品結晶も連続的に取りだす形式である。連続式は大規模な生産に適しているが、結晶の粒度分布を安定させるには結晶の発生や成長速度などを一定に保たなければない。晶析現象は非常に複雑であるため、結晶の発生や成長速度を一定に保つことは難しい。また結晶の発生や成長速度が粒度分布に与える影響は単純ではなく、均一な粒度の結晶を生成するのは困難である。

(10)

一方で回分式は、最初に用意した溶液を装置に供給したら、それ以降は溶液を供給せず、蒸発法や冷却法などによって製品結晶を取り出す形式である。回分式においても、結晶の発生や成長速度が粒度分布に与える影響は単純ではない。しかし、「種晶」と呼ばれる小さな結晶を適切なタイミングで溶液に添加することによって、結晶の発生時期を揃え、均一な粒度の結晶を生成するための方法などの研究が進んでいる（「晶析工学」, pp. 124-131）。製造現場での実務においては、結晶の品質向上、すなわち粒度管理の正確性が求められている。どちらかといえば、連続式よりも回分式で晶析するほうが粒度の管理がしやすい。しかし、生産規模や経済的な事情などから、連続式の晶析でなければならない場合がある。本研究では連続式の晶析装置を対象とする。次節において、本研究で使用している実験装置の仕組みを説明する。

2.2 実験装置

本節では、本研究で使用している実験装置について説明する。はじめに実験装置の構成について説明し、つづいて装置の構成部品について簡単に説明する。装置の構成部品は3 つで、扱う化学物質と溶液、反応槽、およびチューブポンプである。実験装置の構成を説明した後は、この実験装置を動かすための運転プログラムについて説明し、つづいて実際に装置を動かした際に取れる結晶データと、結晶の粒径計測方法について説明する。なお、この実験装置は住友金属鉱山によって作成されたものであり、現在は滋賀大学に設置されている。

2.2.1 実験装置の構成

今回使用する実験装置の構成の全体像は図2.2のようになっている。

(11)

図2.2: 実験装置の構成の全体像まず、装置の一連の動作について説明する。なお、本装置は連続で溶液を供給し続ける連続晶析をおこなう装置である。本実験装置において使用する溶液は、炭酸ナトリウム水溶液と、硫酸亜鉛水溶液の2つである。これら2つの溶液を、純水で洗浄したペットボトルに用意し、チューブポンプを使用して反応槽へと供給する。装置の作動中に溶液濃度を変化させる際は、ペットボトルに刺しているチューブを別濃度の溶液が入ったペットボトルに手動で差し替える。反応槽では炭酸ナトリウム水溶液と硫酸亜鉛水溶液が反応し、炭酸亜鉛の結晶が析出する。また反応槽の底部には撹拌装置があり、溶液をかき混ぜて晶析反応を促している。溶液は反応槽へと連続で供給され続けるため、反応槽の溶液が一定の量となるようにpH・温度計測槽へと適宜排液される。pH・温度計測槽ではpHおよび温度を計測しているのみで、粒径の測定に直接の影響はない。上部の反応槽と同様にして pH・温度計測槽でも、溶液が一定の量となるように適宜排液がおこなわれる。pH・温度計測槽より排液された液体には炭酸亜鉛の結晶が存在するが、これは取り出さない。なぜならば、本実験装置は晶析槽で析出する結晶の粒径を計測するためのものであり、製品を

(12)

製造するためのものではないからである。さいごに、反応槽にて析出する結晶のサンプリングと粒径観測方法について説明する。反応槽で析出した炭酸亜鉛の結晶は、チューブポンプによって溶液とともに汲み上げられ、光学セルを通り再び反応槽へと流れる。光学セル内を流れる炭酸亜鉛の結晶をカメラで撮像し、撮像された画像を画像処理サーバーで処理する。以上が本装置の一連の動作となる。つづいて、扱う化学物質と溶液、反応槽、およびチューブポンプ、これら3つの構成部品について説明する。なお、反応槽内の撹拌装置とチューブポンプにはパルスモーターが取り付けられており、これについてはチューブポンプと共に説明する。

2.2.2 扱う化学物質と溶液の準備方法

扱う化学物質は全部で4つあり、硫酸亜鉛（ZnSO₄）、炭酸ナトリウム（Na₂CO₃）、炭酸亜鉛（ZnCO₃）、および硫酸ナトリウム（Na₂SO₄）である。硫酸亜鉛と炭酸ナトリウムはそれぞれ純水に溶かして硫酸亜鉛水溶液と炭酸ナトリウム水溶液とし、純水で洗浄した2 リットル容器のペットボトルに用意する。炭酸亜鉛と硫酸ナトリウム（水溶液）は、硫酸亜鉛水溶液と炭酸ナトリウム水溶液が反応槽で反応した際に生成される化学物質であり、その化学反応式は式(2.1)である。 ZnSO₄

+

Na₂CO₃

=

ZnCO₃

+

Na₂SO₄ (2.1) 炭酸亜鉛以外はすべて可溶性（水に溶けやすい）であるため、水に溶けてイオンとして存在する。炭酸亜鉛は難溶性（水に溶けにくい）であるため溶液が過飽和状態になりやすく、炭酸イオンと亜鉛イオンが反応すると直ちに溶液中に結晶が析出する。したがって、反応槽で析出する結晶は炭酸亜鉛の結晶のみである。炭酸亜鉛の結晶は白い粉のような結晶である。また硫酸亜鉛水溶液と炭酸ナトリウム水溶液は、それぞれつぎのようにして準備している。硫酸亜鉛水溶液は、0.1mol/lの容量分析用標準硫酸亜鉛水溶液を純水で薄めることで必要な濃度の硫酸亜鉛水溶液を準備している。炭酸ナトリウム水溶液は、286.14g/molの炭酸ナトリウム十水和物1級試薬を純水に溶かすことで必要な濃度の炭酸ナトリウム水溶液を準備している。なお設備の都合上、容積を計測するためのメスシリンダーはなく、電

(13)

子天秤があるのみである。したがって、質量の計測方法はすべてmol/kgとしている。また、準備した2つの水溶液はそれぞれ2リットル容器のペットボトルに入れる。

2.2.3 反応槽

反応槽は図2.3左図のようになっており、反応槽の液面維持と攪拌装置は図2.3右図のようになっている。なお液面維持と撹拌装置は、図2.2にある反応槽とpH・温度計測槽のどちらも同様のものとなっている。図2.3: 左）反応槽,右）液面維持と撹拌装置はじめに図2.3左図について説明する。ペットボトルに準備した硫酸亜鉛水溶液と炭酸ナトリウム水溶液がチューブポンプによって汲み上げられ、反応槽へと連続で供給される。 2つの水溶液が反応槽で反応し、炭酸亜鉛の結晶と、硫酸ナトリウム水溶液（硫酸イオンとナトリウムイオン）が生成される。結晶の粒径を測定するために、反応槽の下部よりチューブポンプによって溶液が汲み上げられ、粒径測定系へと溶液が供給される。粒径測定系に供給された溶液は再び反応槽に返ってきて、反応槽上部より供給される。

(14)

つぎに、反応槽の液面維持と撹拌装置について図2.3右図を用いて説明する。反応槽には、液面の高さを検出するための液面検出電極が設置されている。液面検出電極と、排液に用いられる図2.3右図の左下のチューブポンプは、次のようにして連動している。液面が電極Hに達した場合、左下のチューブポンプが作動して排液が開始される。液面が電極L から離れた場合、左下のチューブポンプが停止して排液が停止される。すなわち、液面は電極Lと電極Hの間に維持されるということである。また、反応槽の底部には撹拌子があり、これはテフロン樹脂で包まれた磁石である。この攪拌子は、反応槽の下にある撹拌用モーターに取り付けられた磁石と連動しており、攪拌用モーターを回転させることで反応槽の溶液を撹拌している。攪拌用モーターの回転速度、すなわち撹拌速度はサーバー上で管理されており、撹拌は事前に決めた一定の速度でおこなわれる。また、攪拌用モーターにはパルスモーターを採用している。パルスモーターについては、次項のチューブポンプとまとめて説明する。

2.2.4 チューブポンプ

本実験装置で溶液の汲み上げや排水をおこなっているポンプは、すべてチューブポンプである。チューブポンプとは図2.4のようなもので、中央にあるローラーがチューブを潰しながら回転することで、チューブ内の溶液を送り出す仕組みになっている。図2.4: チューブポンプ中央のローラーはパルスモーターによって回転させている。パルスモーターとは、ステッピングモーターとも呼ばれるモーターで、パルスという短時間の電気信号を受け取った数の分だけ回転するモーターである。このパルスモーターを使用する利点は、時間当たり

(15)

のパルス数で回転数を正確に制御することができるため、定量性を保つことができる点である。前項の撹拌モーターおよびこのチューブポンプのモーターはともにパルスモーターであるため、撹拌や溶液供給のオンとオフや、撹拌速度、溶液供給速度などを、パルス数を設定することで制御できる。

2.3 運転プログラム

前節では実験装置の構成や、各部品がどのように機能するのかを説明した。本節では、この実験装置をどのようにして動かすのかについて説明する。本実験装置は、装置の傍に設置されているサーバーで動作を管理している。実験装置を動かすにあたり中心となるプログラムは、図2.5のようなものである。図2.5: 運転プログラム

(16)

このプログラムには、溶液供給や撹拌機、粒径計測のためのカメラや照明のオンとオフを、秒単位で記述する。たとえば図2.5のプログラムであれば、「P0撹拌機1」は1秒後に動きだし、133秒後に停止するという指示になっている。このようなプログラムを事前にいくつか作成して組み合わせることによって、常に一定の動作で装置を動かし続けることができる。ただし、溶液の濃度を変化させたい場合はこのプログラムに記述するのではなく、人力での作業が伴う。溶液の濃度を変更するには、溶液の入ったペットボトルに入っているチューブを、手動で別濃度の溶液が入ったペットボトルに移し替えなければならない。だが、この作業は簡単かつ時間がかからないので、実験装置の運転に影響は出ない。注意すべきは、濃度を変化させた時刻はサーバーのデータとして記録されないので、別途記録を取っておく必要があることである。

2.4 粒径の計測方法

本節では、はじめに結晶を撮像するための光学セルやカメラについて説明する。つぎに、撮像された結晶の粒径の計測方法について説明する。さいごに、粒径の計測に関わる注意事項を2つ述べる。図2.6左図にカメラ周りの構図を示し、図2.6右図には実際に撮像された画像を一枚示した。図2.6: 左）カメラ周りの構図,右）撮像画像

(17)

はじめに、光学セルやカメラについて図2.6左図を参照しながら説明する。実験装置を作動させ、光学セルを通過する溶液をカメラで撮像することによって結晶を捉え、画像処理によって結晶の粒径を計測する。光学セルの素材は石英ガラスを使用しており、光の透過や反射による障害が起こりにくいようになっている。また、撮像する瞬間、セルの後方に設置されているLEDライトが点滅することで光量を確保している。撮像された画像はサーバーに送られて画像処理される。つぎに、図2.6右図を参照しながら画像処理による粒径の計測方法について説明する。なお画像処理のロジックは住友金属鉱山が開発したノウハウであるため、詳細な説明は伏せ、大まかな手順についてのみ述べる。セルを通過する溶液を撮像すると、図2.6右図のような画像が得られる。このような画像からピントがあっている結晶近傍の画像を切り取り（図の右上部で赤の四角で囲われている部分）、既定のロジックによって粒径を計測している。さいごに、粒径の計測に関わる注意事項を2つ述べる。 1つは、本研究においては、観測された粒径が50μm以上のデータは除外したことである。なぜならば、撮像された画像に気泡や糸くずなどのゴミが写りこむことがあり、それらを結晶と誤認識して粒径を計測してしまうことがあるためである。粒径が50μm以上のものを気泡やゴミであるとしたのは、試運転（2020年7月頃）より得られた経験から定めた値である。もう1つは、実験装置の作動中に光学セルが曇ってしまうことである。実験装置を連続で作動し続けていると、光学セルに結晶がこびりついて、セルが曇っていくのを目視で確認することができる。セルが曇っていくと図2.6右図のように撮像される画像の輝度が低下していく。セルの曇りが激しい場合は、クエン酸や酢酸などで洗浄することで曇りを解消することができる。

(18)

第

₃

章

実験条件と実験結果

本章では、本研究の分析において使用する実験データを取得した際の実験条件と、その実験結果について説明する。

3.1 実験条件

本節では、本研究で分析するデータを取得したときの実験条件について説明する。

3.1.1 実験条件とは何か

はじめに、実験装置の運転条件と実験条件の区別について説明する。運転条件とは、実験装置を作動させるときに操作や監視をする変数である。たとえば、溶液の攪拌速度や溶液の供給速度などのことである。一方で実験条件とは、何らかの目的のために、あらかじめ運転条件をある特定のものにすると決めた条件のことである。たとえば、実験装置の作動中に、運転条件の1つである溶液濃度を濃くしたときの粒度の変化を見るという目的があった場合、「濃度を1から2に上げる」というようなものが実験条件である。なお、本実験装置で考慮される運転条件は複数あるが、主に考慮すべきものは5つあり、温度、pH、反応槽の溶液攪拌速度、硫酸亜鉛水溶液と炭酸ナトリウム水溶液の溶液供給速度およびそれらの溶液濃度である。

(19)

3.1.2 本研究における実験条件

はじめに実験条件を決めるにあたり考慮すべきことについて説明した後に、本研究で定めた実験条件について述べる。本研究では、運転条件を変化させたときに粒度がどのように変わるのかを分析することにより、結晶品質の安定化方法を探るということが目的であった。変化させることができる運転条件は複数あるが、本研究では硫酸亜鉛水溶液と炭酸ナトリウム水溶液の溶液濃度を変化させることにした。溶液濃度を変化させる理由は、複数ある運転条件のうち、溶液濃度を変化させるのが容易であるから、かつ最も粒度の変化が大きいと考えられるからである。なお、溶液濃度の変化による影響のみを調べるために、その他の運転条件である温度や溶液供給速度、攪拌速度などはすべて一定にしている。溶液濃度の値や変化させる範囲については、事前におこなった試運転（2020年7月頃）の結果より、硫酸亜鉛水溶液と炭酸ナトリウム水溶液ともに0.001mol/kg前後から0.004mol/kg前後の間の濃度が適切であると判断した。濃度範囲の判断は、下限濃度については分析に必要な結晶の観測数から考慮し、上限濃度については実験装置の画像処理できる限界を考慮したものである。また同様にして試運転の結果より、以下のことを決めた。実験装置の作動中にセルの洗浄をおこなうと、観測されるデータに大きな乱れが発生することが確認されたので、実験中にセルの洗浄はおこなわないこととした。上述したことを踏まえて、実験条件を以下のように設定した。（実験条件）硫酸亜鉛水溶液と炭酸ナトリウム水溶液の濃度を、ともに0.001mol/kgから 0.0015mol/kgへ変化させるなお、溶液濃度を変化させるのは、観測される粒度がある程度定常状態に落ち着いていると見られたときである。定常状態に落ち着いているかどうかは、装置の運転中にサーバーからリアルタイムで粒径のデータを抽出し、粒径から算出したD50の変動が小さくなったかどうかで判断した。

3.2 実験結果

前節で定めた実験条件に従って実験装置を動かし得られた実験結果を、図3.1に示す。

(20)

図3.1: 実験結果

はじめに図3.1の見方について説明し、次にこの実験結果から読み取れることについて述べる。

この実験は2020年9月4日におこなったものである。変数は上から順に、画像平均輝

度（avg_kido）、観測結晶数（count）、そしてD50である。横軸の lt_timeは観測された

時間であり、プロット上のひとつひとつの点はロットと呼ばれる1分ごとの統計値を示している。言い換えれば、プロット上のひとつひとつの点において、図2.1のような粒度分布が観測されており、その統計値のみを可視化している。16:45にある黒の縦線は濃度を変化させた時刻を示しており、硫酸亜鉛水溶液と炭酸ナトリウム水溶液の濃度をともに 0.001mol/kgから0.0015mol/kgへと変化させている。図3.1から読み取れることは複数ある。たとえば、運転開始時はD50が安定するまで時間がかかることや、画像平均輝度が低下し続けていることが読み取れる。これらの現象は、前の運転時の反応槽の状態を引きずりD50がしばらく安定しないからであったり、光学セルに結晶がこびりつくこと（目視によって確認ができる）で輝度が低下しているからだと考えられる。図3.1から読み取れることで特に重要なのは次の2つのことで、（１）濃度を濃くした時

(21)

点より観測結晶数が増加し、D50も大きくなっているということと、（２）輝度の低下、あるいは観測結晶数の低下がD50にどのような影響を与えているかがわからないことである。なぜなら、運転条件である溶液濃度を変化させた際にD50が変化したということは、粒度分布に対して何らかの影響があったということである。また、輝度や観測結晶数の低下とD50の関連性について考慮しなければ、溶液濃度の変化がD50に与えた影響のみを調べることは困難だからである。次章では、この実験結果を用いて濃度変化の粒度への影響を調べるために、状態空間モデルを用いた分析をおこなう。

(22)

第

₄

章

状態空間モデルによる分析

本章では、3章で得られた実験データより、運転条件である溶液濃度を変化させたときに粒度分布がどのように変わるかを、状態空間モデルを用いて分析する。はじめに、濃度変化や輝度低下、観測結晶数の低下などが粒度分布にどのような影響を及ぼしているのかを、探索的データ解析によって調べる。つづいて、探索的データ解析より示唆されることより仮説を立てて、状態空間モデルを構築する。構築したモデルのパラメータ推定には MCMC法を用いる。

4.1 探索的データ解析

本節では、溶液濃度を変化させたときの粒度分布の変化を調べるために探索的データ解析をおこなう。探索的データ解析とは、アメリカの統計学者ジョン・テューキーによって提唱されたものである[7]。探索的データ解析において重要な観点は、次のようなものである。すなわち、置かれた前提から導出される命題群が形づくる演繹的体系も大事であるが、現実のデータを注意深く観察し、数学的理論体系をよりどころとしたデータに基づく推論を実行する、という観点である[8]。本研究においては、晶析工学の知見から演繹的に導き出された数理モデルではなく、実験によって得られたデータに基づく推論から数理モデルを構築している。したがって、本節でおこなっている探索的データ解析は、本研究において非常に重要な部分である。

(23)

4.1.1 粒径域で分割した粒度分布

3章の最後に、図3.1から読み取れる重要なこととして次の2つのことを挙げた。（１）濃度を濃くした時点より観測結晶数が増加し、D50も大きくなっているということと、（２）輝度の低下、あるいは観測結晶数の低下がD50にどのような影響を与えているかがわからないことである。（１）に関しては、次の2つのことがいえる。1つは、濃度を濃くした時点より観測結晶数が増加したことについてである。晶析の理論的にも、溶液濃度を濃くすると、溶液に含まれる炭酸イオンと亜鉛イオンの量が増加するため結晶数が増えるというのは確からしいことである。もう1つは、D50が大きくなっていることについてである。D50が大きくなるということは、結晶群が全体的に大きくなっているか、あるいは大きな結晶が増加しているかであると考えられる。なぜD50が大きくなっているのかは図3.1のプロットから読み取ることはできないため、データの見方を変える必要がある。また（２）の輝度の低下や観測結晶数の低下がD50へ与える影響に関しても、図3.1のプロットからは読み取ることはできない。溶液濃度を濃く変化させたとき、D50が大きくなっているのはなぜかということ、さらに輝度低下や観測結晶数の低下とD50（粒度分布）の関係性を確認するために、図4.1のようにして粒度分布を11の粒径域に分割した。粒度分布をこのように分割したのは、溶液濃度を濃く変化させたとき、粒径の大きさによって変化が異なるかどうかを確認するためである。したがって、分割した基準となる粒径に特別な意味はない。なお粒径域を分割する基準となる粒径は、粒度分布の10%ごとの分位点となる粒径をおおよその目安に定めたもので、各粒径域の大きさの範囲は表4.1にまとめた。

(24)

表4.1: 粒径域とその大きさ粒径域（grp）下限上限 u4.5 0μm 4.5μm u4.75 4.5μm 4.75μm u5 4.75μm 5μm u5.5 5μm 5.5μm u6 5.5μm 6μm u6.5 6μm 6.5μm u7 6.5μm 7μm u8.5 7μm 8.5μm u10 8.5μm 10μm u13 10μm 13μm u50 13μm 50μm 図4.1: 粒径域で分割した粒度分布

(25)

実験によって得られたデータで、分割した粒度分布のそれぞれの粒径域に属する結晶数の割合をロットごとに集計し、その割合の遷移を示したものが図4.2の下部のプロットである。図4.2: それぞれの粒径域に属する結晶の割合の遷移下部のプロットの縦軸のrateは、それぞれの粒径域に属する観測結晶数の割合を示している。ロットでの集計時点は上部のプロットと同じであるが、割合の遷移は折れ線グラフのみで表示している。図4.2からは、先ほど図3.1から読み取ることができないと述べた2点それぞれに関して言及できることがある。

(4.2) 式(4.1)と式(4.2)はともに方程式の形で表現されることが多く、状態の変化を表す方程式は「状態方程式」、観測値を表現する方程式を「観測方程式」と呼ぶ。また、本研究において状態空間モデルを用いたのは、次の理由からである。すなわち、分割した粒度分布の各粒径域の割合について、直接観測できない反応槽の真の割合が濃度に依存して遷移していく様子を、確率的に観測される各粒径域の割合から推定するためである。

(28)

4.2.1 粒径域ごとに溶液濃度に依存した結晶成長の状態空間モデル

前節でおこなった探索的データ解析より、次の2つの考察が得られている。1つは、観測結晶数が低下していく状況においても、粒度分布を適当な粒径域に分割した場合の、粒径域ごとの割合は安定していること。もう1つは、濃度を濃く変化させたときに大きな結晶の割合が増加したのは、ある粒径域から次の粒径域に結晶が成長するスピードに変化が起こったからであると推察されることである。本研究において構築した状態空間モデルは、上記2つの考察を加味し、分割した粒径域ごとに、溶液濃度に依存して結晶成長割合が異なるというモデルである。はじめに、状態モデルと観測モデルおよびモデルのイメージ図と、以降で用いる変数と記号の意味をまとめる。つづいて、モデルの構造について（１）状態モデル、（２）観測モデルの順に説明する。 x_𝑡+1

= 𝐴(𝑐

_𝑡

)

x_𝑡 (4.3) y_𝑡

∼ 𝑀 (

x_𝑡

, 𝑛

_𝑡

)

(4.4)

𝑝

_𝑡,𝑖

_𝑡,1

, … , 𝑦

_𝑡,11

)

𝑡：

𝑦

_𝑡,𝑖 は

𝑡

𝑖

の観測結晶数

𝑀

：多項分布

𝑛

_𝑡：

𝑡

時点における観測結晶数

𝑝

_𝑡,𝑖：

𝑡

𝑖

の結晶成長割合

𝑞

：放流率

の割合で放流され、放流されたぶんは一番小さな粒径域 1に加えるというは次のような理由からである。放流によって減少した割合がすべて粒径域1へ移動するとしたのは、溶液より析出する結晶は、最初から大きな結晶として現れるのではなく、はじめは小さな結晶として現れると考えられるためである。なお、粒径域の割合の総和が必ず1となるように、遷移行列

𝐴

の各列で総和が1となるようにしている。ここで成長割合

𝑝

_𝑡,𝑖 について説明する。

𝑝

_𝑡,𝑖 は、

𝑡

4.3.2 MCMC

法による計算

MCMCとはMarkov Chain Monte Carlo（マルコフ連鎖モンテカルロ）の略称で、これ

は乱数発生のためのアルゴリズムである。この乱数発生アルゴリズムであるMCMCを事後分布の計算に応用したものを、ここではMCMC法と呼ぶことにする。なぜ事後分布の計算にMCMCを応用するかは、次のような理由からである。ベイズの定理より、事後分布の確率密度関数は、尤度関数と事前分布の確率密度関数の積を、尤度関数と事前分布の確率密度関数の積を積分したもので除したものである。しかし、この事後分布の確率密度関数は、モデルが複雑となることが多く、特にパラメータが多い時には非常に複雑になり、積分計算が困難になる。そこで、MCMCによって事後分布の確率密度関数に従う乱数を発生させることで、パラメータがある範囲の値を取る確率や、事後分布の期待値などが計算できるようになる。 MCMC法による計算方法には、MH法（メトロポリス・ヘイスティング法）やHMC法

(33)

（ハミルトニアン・モンテカルロ法）などがある。MH法とは「仮の数値をそれぞれの変数に投入し、他の変数値を固定したときにある変数の2つの数値を比較して、データの当てはまり（尤度、もっともらしさを表す統計値）と事前分布の確率の積が高い方を採用していく手順を取って、変数値の分布（事後分布）を描き出す方法である（計算の過程で、ある一定の割合で尤度と事前確率の積が低い値も採用するようになっている）」[12]。HMC 法の基本的な考え方はMH法と同じであり、Stanで使用するのはHMC法である。 MCMC法による計算の結果、抽出される値がある事後分布に収束しているかどうかを判断するために、収束診断が必要である。この収束診断法方法は複数あるが、本研究においては次の2つの方法で収束しているかどうかの判断をおこなった。1つは、トレースプロットと呼ばれる時系列プロットの目視による判定である。MCMC法では初期値の異なる複数の連鎖（chains）で乱数を発生させることができ、これを時系列でプロットしたものがトレースプロットである。トレースプロットですべての連鎖が別れることなく交わっていることが収束の判断基準である。もう1つはRhatが1.1未満になっているかどうかである。Rhatは、事後分布からのサンプリングが収束したかを見る指標で、1へ収束する指標である。初期値の異なる連鎖の収束先が異なっていた場合などは、Rhatは大きくなる。すべての

𝑘

₁と

𝑘

₂、

𝑞

、および初期値x₁について、Rhatが1.1未満であることが収束の判断基準である。

4.3.3 Stan

の準備

MCMC法を実行するためのソフトウェア環境にはBUGSやStanなどがあるが、本研究ではStanを利用する。Stanに受け渡すためのデータは、すべてRで準備している。R で生データが保存されているサーバーのデータベースから必要なデータを抽出し、モデルの推定に必要なデータセットを作成した。StanはRStan[13]というパッケージを用いて、 Rから呼び出している。 MCMC法の計算にあたり、パラメータそれぞれに次のような制約を設けた。放流率

𝑞

に関しては、排液によって反応槽から結晶がなくなる割合なので下限を0とし、上限を 0.4とした。上限を0.4としたのは、単位時間あたりに、反応槽にある結晶がほぼ半分なくなることはないと考えられるためである。ミカエリス・メンテン式のパラメータ

𝑘

₁ は、溶液濃度なので下限を0とし、上限は1とした。3章で述べたように、本実験装置においては0.004mol/kg程度が実験に適した濃度の上限であるため、上限を1としても推定に支障はないと考えられる。ミカエリス・メンテン式のパラメータ

𝑘

₂は、成長割合（の上限）

(34)

なので下限を0とし、上限は1とした。つぎに、MCMC法を実行する際の計算回数などの指定について説明する。MCMC法では、アルゴリズムの都合上、乱数が初期値へ依存してしまうことや、連続する採用値に自己相関が含まれてしまうことが起こりえる。こういった問題を回避するために、次のような対策をする。はじめにウォームアップ期間（warmup）を設定する。ウォームアップ期間とは、生成された乱数を切り捨てる期間である。ウォームアップ期間を指定することで、初期値への依存性を下げることができる。つぎに間引き間隔（thin）の設定をする。間引き間隔とは、生成された乱数を間引く間隔の指定である。たとえば間引き区間を2とすると、生成された乱数のうち1つ目の値を採用し、2つ目の値を切り捨て、3つ目の値を採用するといったようになる。乱数を間引いて採用することで、乱数の自己相関が起きることを防ぐ。さいごにマルコフ連鎖の数（chains）を設定する。マルコフ連鎖の数とは、初期値の異なる乱数の連鎖を発生させる数の指定である。異なる初期値から出発した乱数の連鎖が交わっているかどうかを確認することによって、事後分布が収束したかどうかの目安となる。本研究におけるMCMC法の計算は、20万回おこない、初期値からの10万回を切り捨て（warmup=100,000）、以後の100回ごとに結果を抽出（thin=100）する計算を初期値の異なる3本の連鎖（chains=3）についておこない、事後分布を求めた。すなわち結果として得られる事後分布からの標本の大きさは3千である。なお、サンプリング結果が再現できるよう乱数のシードを1229とした。

(35)

第

₅

章

分析結果と考察

本章では、はじめにMCMC法による計算結果の収束診断をおこなう。つづいて、推定されたパラメータを用いて、4章にて構築した状態空間モデルによって生成されるデータと観測データを比較し、モデルの妥当性を確認する。

5.1 MCMC

第

₆

章

今後の課題と結論

6.1 今後の課題

本節では、本研究において未検証となってしまったことや明らかにしきれなかったことなどについて述べる。今後の課題として挙げられるのは、次の7つの事項である。 1. 本研究におけるモデルの検証 2. 適切な粒径域の分け方 3. 輝度低下が起きない実験装置の開発 4. ミカエリス・メンテン式のパラメータ推定 5. 状態モデルの過程誤差の組み入れ 6. 初期値の推定 7. 晶析メカニズムの解明１．本研究におけるモデルの検証本研究で構築した状態空間モデルは、2020年9月4日におこなった実験のデータから考案し、そのデータに対してのみパラメータ推定をおこないモデルの評価をした。このモデルが妥当であるかどうかは、次のようにして検証する必要があるだろう。はじめに、 2020年9月4日におこなった実験と同じ実験条件で再実験をおこない新たなデータを計測する。つぎに、本研究で構築したモデルおよび推定したパラメータを用いて、再実験より得られた新たなデータを再現できるかどうかを確かめる。しかし、晶析現象の複雑さから、実験条件をすべて同じにして再実験をおこなったとしても、2020年9月4日と同じ

連続晶析における粒度分布遷移過程の状態空間モデルによる分析

＜ 修 士 論 文 ＞