Microsoft Excelソルバーによる線形計画法の滴定曲線解析への応用

吉村 季織, 岡崎 正規, 中川 直哉


Return

1 緒言

パーソナルコンピュータ上で用いられている表計算ソフトウェアの中でも、Microsoft Excel(以下Excel)は現在最も広く利用されている。Excelには通常の表計算ソフトとしての機能以外にも、統計ツールやデータベース機能など様々な機能が付加されている。ソルバーもそのひとつであり、専門的な知識を持たなくても最適化を行うことが可能となる。
カーブフィッティングを行うとき、しばしば最小二乗法(LS)が用いられている。前報では、ソルバーによる非線形最小二乗法による滴定曲線の解析、pK及び濃度の算出を行った[1]。LSでは、誤差の二乗和を最小化させるが、誤差の絶対値の和を最小化させることで、最もデータに近づいた予測曲線を見出すことができる。予測曲線の理論式が線形で表されるなら、このフィッティングの問題は線形計画法(LP)を用いて解くことができる。
本稿では、まずLPによるカーブフィッティングを、Excelソルバーを用いて解く方法について解説する。さらに応用として、最近の研究でも用いられている、pK分布[2 - 5]を作成する。

2 LPによるカーブフィッティング

LPとは、目的関数、制約式がともに線形で表される場合である。一般的な場合を示すと、変数viに関して、

となる。ajibicjは係数。最大化問題は、目的関数の右辺の符号と、制約式の不等号を逆転させることで最小化問題となる。LPで解くときにまず考えなくてはならないのは、いかにして解きたい問題を線形で表すかである。
m個の独立変数(x1, x2, ..., xm)からなる関数y' = S(cixi)の場合を考える。ciは係数である。j番目の実際のデータ(xj1, xj2, ・・・, xjm, yj)と関数y'による予測値yj'の誤差の絶対値の和S|yj' - yj|を最小化すれば、カーブフィッティングができる。しかしながら、この和は線形ではない。ここで、誤差の大きさをvjとし、変数として扱うことにより、上記問題は以下のように書き換えられる。

ここで、vjは実測値と予測値の誤差を直接求めているのではなく、あくまで変数であり、制約式によって誤差の大きさになるよう調整されている。目的関数fは、明らかに線形である。関数y'は線形なので、すべての制約式もまた線形である。ゆえにこの問題はLPによって解くことができる。この場合、関数y'に含まれるすべてのcivjが変数となる。目的関数の中にはciが無いが、これはciの係数がすべて0であると考えるとよい。この問題を、たとえばC言語の関数で解くには、その関数で解くことのできる入力形式や、式(1)のように一般化された形にする必要がある。特に行列の生成は面倒である。しかしながら、Excelソルバーを用いれば行列のことを考えることなく、ほとんど式(2)の状態のままで問題解決まで進めることができる。
最もよく用いられる線形回帰へのLPの適用を例としてあげる。まず、Figure 1の左上に示すようなワークシートを用意する。セルB5からC14はこれからフィッティングするためのデータである。F1およびF2には、回帰直線の係数abとし、初期値として0を与えた。E5からE14はyの予測値y'で、たとえばE5には"=$F$1*B5+$F$2"が入力されている。F5からF14はデータとの誤差vで、式ではなく実際に値を入れなくてはならないが、LPでは変数として扱われるので、初期値として0に設定した。また、F5からF14間での合計値(=SUM(F5:F14))をF15に入力した。これが、目的関数に相当する。G列、H列はそれぞれ制約式に当たる。5行目を例に挙げると、セル内の記述内容は、それぞれ"=E5+F5"、"=E5-F5"である.


Figure 1. Preparation to solve a linear regression problem by LP in Excel solver.

次に、回帰のためのソルバーの設定を行う(Figure 1 下段)。ソルバーの設定とワークシートの関係を、矢印で示した。最小化の目的となるのは、Svだから、目的セルには"$F$15"を入力する。最小化問題であるから、目標値は"最小値"とした。変数は、回帰式の係数とすべてのvなので" $F$1:$F$2,$F$5:$F$14"を入力した。10個のデータがあるので、30個の制約条件が存在するが、同一内容の条件についてはFigure 1に示したように、"$G$5:$G$14 >= $C$5:$C$14"のようにまとめて設定することができる。ここで、実行をすることで最適化は可能であるが、ソルバーには線形計画法を解くためのモジュールが搭載されているので、今回の場合それを利用した。パラメータ設定で、“オプション”を押すと、Figure 1右図のようなオプション設定がでてくる。この中の、“線形モデルで計算”をチェックした。その他の設定は変更しなかった。これで設定は終了で、後は実行すればよい。


Figure 2. Comparison of results from linear regression by LP

実行した結果、回帰係数abそれぞれ、1.8、5.6が得られた。Figure 2にデータと得られた回帰直線(実線)の関係を示した。比較として最小二乗法による回帰結果も示した。これらの2本の直線は、わずかに異なっているが、よく似た結果となっている。データのばらつきが大きくなるにつれて、2本の差も大きくなってくる。しかし、回帰の方法が異なっているので、どちらがより勝っている方法かを示すものではない。

3 LPによる滴定曲線解析

3. 1 滴定曲線の線形変換

前述したようにLPで問題を解くためにまず必要なことは、問題を線形で表すことである。一般に用いられる滴定曲線は(滴定体積)-(pH)で表されるが、これは非線形である。そこで、LPで扱えるようにするために、解離酸総濃度(total concentration of dissociated acids : TCDA)[6]を導入し、以下の式にしたがって、滴定曲線を(pH)-(TCDA)形式に変換する。

ここで、Vaは試料溶液の初期体積、cbVbはそれぞれ滴定剤の濃度と、滴下体積である。(3)式で求められるTCDAは実験で得られる値のみを用いるので、TCDAexpと表す。
次に、多価酸は異なる多数の一価酸の等濃度混合物として見なされ、それぞれの一価酸の解離定数は多価酸の逐次解離定数と等しくなる、という仮定からTCDAを理論的に導く(TCDAcal)と、

となる。ここで、ciはpKiを持つ酸の濃度である。式(4)はpHに関しては非線形であるが、pKiをあらかじめ与えることによって、Xiに関して線形と見なすことができる。すなわち、あらかじめ適当なpKiを与えて、TCDAexpy、TCDAcaly'に置き換えてやれば、式(3)と(4)は式(2)で表されるLPにすることができる。このLPを解くことで最適なciの組を見つけることができる。
しかしながら、未知試料中に含まれている酸の種類数と、それぞれの酸のpKが一般に不明である。このようなときは、適当な範囲のpKを適当な間隔で区切ってあたえてやるとよい[6]。一般に(強酸、強塩基以外の)酸のpKは、0から14の間に存在する。また、一般に得られる滴定曲線もこのpKの範囲を外れた酸の正確な情報を含んでいない。そのため、pKの範囲を0から14とした。また間隔を0.2とした。すなわちpK1 = 0.0, pK2 = 0.2, ..., pK71 = 14.0である。

3. 2 滴定実験

5.12 mMのクエン酸水溶液100mlを135 mMの水酸化ナトリウム水溶液で滴定した。滴定体積は1滴あたり0.10mlとし、pH10.03まで行い116組の実験値を得た。


Figure 3. Original titration curve(upper) and converted titration curve(lower) for 5.12 mM citric acid.

Figure 3に滴定曲線と(pH)-(TCDA)型に変換した滴定曲線を示した。

3. 3 LPのための準備

Figure 4に滴定曲線を解析するためのワークシートの一部を示す。A列及びB列の5行目以降が実験値で、その値をもとにD列でTCDAexpを算出した。2行目F列からBX列までは設定したそれぞれのpKを入力し、それぞれの列に対応する3行目がその酸の濃度とした。濃度は、LPにおいて変数となるので、初期値を0とした。それぞれの列に対応する5行目以降で、実験によって得られた各pHにおけるciXiを算出した。絶対参照の使い方に注意すると、簡単に全体へコピーできる。それぞれのpHに関して、SciXiすなわちTCDAcalをBY列5行目以降で計算した。BZ列にTCDAexpとTCDAcalの差の絶対値vとし、CA列、CB列でそれぞれTCDAcal + v、TCDAcal - vを計算した。最後に、BZ列3行目でvの総和を計算した。


Figure 4. Preparation of work sheet to determine the concentrations and pK values of citric acid.

ソルバーの“パラメータ設定”及び“オプション設定”は、前述の例に少し変更を加えた。前述の例では、変数は回帰係数abと誤差の絶対値viであり、viのみに非負制約を設定した。滴定曲線解析の場合、変数は回帰係数に当たる濃度ciと誤差の絶対値viである。濃度は当然非負数のみをとるので、今回の例ではすべての変数に非負制約を設定する必要がある。そこで“パラメータ設定”の“制約条件”で設定せずに、“オプション設定”の“非負数を仮定する”をオンにした。これによってすべての変数が非負として扱われる。ここで注意しなくてはならないことは、Excelのソルバーで扱うことのできる変数の数が200までに制限されていることである。滴定曲線の解析では、滴定データの数と設定したpKの数の和が変数の数となる。本例では、滴定データ116点、71pKなので、合計187変数でありソルバーで扱うことができる。もし200を越えるようならば、変数を減らす工夫をする必要がある。さらに“オプション設定”において反復回数を1000に増やした。これは、初期値の100では最適化が終了せずに計算途中で継続するか否か確認させられるので、これを回避するためで、計算精度には影響しない。
計算は、Windows NT 4.0上で動作している、Excel 2000にて行った。

4 結果と考察

計算により得られた数値結果を折れ線グラフとして、Figure 5に示した。このような図をpK分布、pKスペクトルなどと呼ぶ。Figure 5より、pK3.0、4.5、5.5付近にそれぞれ大きなピークが現れ、pK7付近と、10と11の間に小さなピークが見られた。前者の3本のピークがクエン酸由来で、それぞれクエン酸の逐次解離に対応する。後者の小さなピークは、炭酸などの不純物により現れたものと思われる。同じ滴定曲線をシンプレックス法を用いて作成した自作プログラム[7]で解析し、ソルバーと比較した。シンプレックス法は非負制約付線形計画法を解く最も一般的な方法である。前者3本のクエン酸由来のピークでは、ソルバー算出値と自作プログラム算出値間での相対誤差の平均が、pKでは3.10×10-8 %、濃度では8.71×10-7 %と非常に小さく有意なものではなかった。pK7付近ピークでも0.02 %以下の誤差であり、有意とは言えない。pK10-11のピークは、自作プログラムでは1本しか現れなかった。これは計算法によってpKが変化する不安定なピークであることから、比較ができなかった。このように、Excel ソルバーで解いた結果は、専用のプログラムで解いた結果と遜色がなかった。
それぞれのピークは、1つ以上のpKとそれに対応する濃度cによって構成されている。例えば、pK3.0付近のピークは、pK2.8の酸3.16 mMとpK3.0の酸2.09 mMからなる。pK5.5付近のピークが大きく見えるのは、このピークに含まれるpKが1つだけであるためである。また、ピークとして認められるpK以外の領域では、すべてのpKiに対応する濃度ciは0.0 mMであった。TCDAexpとTCDAcalの差の平均は、1.37×10-2 mMと非常に小さく、良くフィッティングが行われていることが分かる。


Figure 5. pK distribution of 5.12 mM citric acid solution

pKの間隔が0.2と広いため、それぞれのピークも幅が広くなった。しかしながら、1種類の解離に対応するpKは理想的には1つであり、そのことをpK分布で表すと、完全な1本線になるはずである。ピークの出た領域内をさらに細かいpK間隔で区切り、同様の計算を繰り返すことで、ピークはより細くなり、より正確な濃度やpKを見積もることができるようになっていくことが明らかになっている[6]。このことは、ピークの現れている領域内に真のpKがあることを示唆している。今回得られた結果は荒いものであるが、0.0でないciに対応するpKiを加重平均により補正することでより正確なpKを見積もることができる。

濃度は、ピークを構成するpKに対応するcの合計とする。

これらの補正法によって得られた、クエン酸のpKと濃度をTable 1に示した。補正によって算出されたpKの値は、文献値に良く一致していることが分かる。調製した濃度は5.12 mMであったから、やや誤差が大きいようである。この誤差は試料についての情報が全くなかったことや多価酸を1価酸の混合で近似したことよりも、pHメーターの持つ誤差など、実験時の誤差の方が大きく影響している。しかしながら、この3つの酸の平均濃度は5.10 mMと調製した濃度に充分近く、等量点の見積もりは正確に行うことができる。

Table 1. pK value and concentration of citric acid.
pKconcentration
mM
reference[12]experimental
2.902.885.25
4.344.345.18
5.665.604.87

5 終わりに

Excelソルバーは、Excelユーザーなら誰でも使うことのできる最適化ツールである。本稿では、ソルバーに内包されている線形計画法のモジュールを活用して、水溶液未知試料の滴定曲線より、pK分布を見積もる方法を示した。もしこの手順を自作プログラムとして作成するならば、線形計画法の解を求めるモジュールを作らなくてはならないほか、線形計画法に則るために目的関数の係数ベクトル、制約式の係数行列などを生成するためのルーチンを記述しなくてはならない。もちろんソルバーで問題を解くとなると、計算時間を犠牲にしなくてはならないが、滴定に要する時間に比べると充分に小さく、実験の中での律速段階にはならない。
本法は、最低限の試料の情報のみで、pK分布を見積もるので、特定の酸の定量やpK決定など、精密な分析を行うことには向いていない。逆に、pK分布の傾向を見積もることができるため、未知試料や複雑な混合系の分析に適している。実際に、pK分布は酸化物[8]や粘土鉱物表面[9, 10]やイオン交換樹脂[11]、微生物表面[4, 5]、有機物[3]など、不均一表面、複雑な系の分析に使われている。
最近のパソコンには、あらかじめExcelがインストールされていることが多いため、本法はあらゆるユーザーに、pK分布を見積もる機会を提供するものである。

参考文献

[ 1] 吉村季織, 岡崎正規, 中川直哉, J. Chem. Software, 7, 191 (2001).
[ 2] P. Brassard, J. R. Kramer and P. V. Collins, Environ. Sci. Technol., 24, 195 (1990).
[ 3] D. S. Smith and J. R. Kramer, Environ. Int., 25, 307 (1999).
[ 4] J. S. Cox, D. S. Smith, L. A. Warren and F. G. Ferris, Environ. Sci. Technol., 33, 4514 (1999).
[ 5] I. Sokolov, D. S. Smith, G. S. Henderson, Y. A. Gorby and F. G. Ferris, Environ. Sci. Technol., 35, 341 (2001).
[ 6] N. Yoshimura, M. Okazaki and N. Nakagawa, Anal. Sci., 16, 1331 (2000).
[ 7] 吉村季織 未発表
[ 8] C. Contescu, J. Jagiello, and J. A. Schwarz, Langmuir, 9, 1754 (1993).
[ 9] T. J. Bandosz, J. Jagiello and J. A. Schwarz, J. Colloid Interface Sci., 182, 570 (1996).
[10] T. J. Bandosz, J. Jagiello and J. A. Schwarz, J. Phys. Chem., 100, 15569 (1996).
[11] S. Y. Bratskaya and A. P. Golikov, J. Anal. Chem., 53, 234 (1998).
[12] 日本化学会編, 化学便覧 基礎編II, 改訂4版, 丸善 (1993), p.318.


Return