回帰計算

ExcelのソルバーやPythonの最適化を使った回帰計算例、回帰計算の注意点等を整理していきます。

NRTLのbinary parameterのRegression例

 

違いを見やすくするために、現在のbinary parameterがほぼ理想状態からスタートし、水-エタノールの回帰を行います。

 before

Excelのソルバーを使った場合

 

実験データと計算上の値との自乗誤差を積算し、誤差が最小となるパラメーターをExcelのソルバーで求めます。

 afterExcel  

以下はwater-methanolの例になります。

 

(NRTLの詳細については、こちらを参照ください)

Pythonの最適化関数を使った場合

 

scipy.optimizeのleastsqを使っていきます。

回帰計算時の注意点

Anscombe‘s Quartet(アンスコムの例)

 Anscomb

・。

シミュレーション前のデータ処理

Data Reconciliation

最小二乗法と重みづけ

プラントの生データは、流量計の誤差や成分分析の測定誤差等から、ままマテバラが取れていないことがあります。

プロセスシミュレーターで計算を行う前に、マテバラが取れたデータを用意する必要があります。

以下Excelを使う方法とPythonを使う方法をご紹介します。

データ補完

Cubic Spline

お勧め書籍:Cによる理工学問題の解法、佐藤 次男 (著), 中村 理一郎 (著), 伊藤 惇 (著)

Natural Cubic Splineの考え方

n個のデータ(x1,y1)~(xn,yn)について、以下3つの条件を満たす3次スプライン関数S(x)を考えます。(x1 < x2 < ... < xi < ... < xn)

①1次導関数S'(x)が連続であること : \(S'_i(x_{i-1})=S'_i(x_{i})\)

②2次導関数S''(x)が連続であること : \(S''_i(x_{i-1})=S''_i(x_{i})\)

③i=1~nに対し、S(xi)=yiであること

さらにNatural Cubic Splineの場合、境界条件として、S''(x)の始点と終点を0としますが、こちらは後述します。

まず条件②から考えていきます。(大概の解説はS(x)=\(ax^3+bx^2+cx+d\)から追いかけていきますが、私は上記お勧め書籍で説明されている2次導関数から追いかける方が分かりやすいので、②からスプライン関数を求めます。)

3次関数S(x)の2次導関数は、1次関数となります。

ここで、S''(x)=M(x)、i=1~nに対して、Mi=M(xi)とすると、前述の2次微分の連続性から、\(x_{i-1} \leqq x \leqq x_i\)において、下図のようにS''(x)を表現することができます。

 Interpolation1

S''(x)を積分し、S'(x)、S(x)を求めます。

 Interpolation2

S(x)が求まったので、続いて条件③を評価します。\(S(_{i-1})=y_{i-1}\)、\(S(_{i})=y_{i}\)より、以下積分定数C1、C2を求めます。

 Interpolation3

最後に、条件①を評価します。

S'(x)が連続であることから、\(x_i\)において、\(S'_i(x_i)\)と\(S'_{i+1}(x_i)\)は等しくなります。これを整理すると、

 Interpolation4

上式に前述のNatural Cubic Splineの始点と終点の境界条件である\(M_1=0, M_n=0\)を加え、\(M_1 ~ M_n\)が求まります。

上式の\(M_2 ~ M_{n-1}\)を求めるには、三重対角行列アルゴリズムであるトーマスアルゴリズム(TDMA)を用います。

以下具体的な解き方については、「Silverlightでツール作成」の項で後述いたしますが、Excelを使う場合は、逆行列が簡単に求まりますので、もっと簡単です。

Scipyのinterpolateを使う方法

scipyのモジュール:interpolate.CubicSplineを使うことでcubic splineの作成が可能です。

Excelのminverse、mmlutを用いる場合

例として、次の6点を通るスプライン曲線を考えます。

(11.1,34.3),(20.9,62.1),(40.0,64.5),(65.6,75.9),(73.9,90.7),(90.9,95.9)

x1~x6、y1~y6を代入し、M2~M5を求めていきます。(natural cubic splineの条件から、M1=M6=0)

 Interpolation5

Excelの場合は、minverse関数で逆行列を、mmult関数で行列の積を簡単に求めることができます。

 Interpolation6

試しに元データの中間点のXの値からYの値を求めプロットすることで、以下のようにCubic Splineが求まります。(もっと刻みは増やした方が良いです)

 Interpolation7

Silverlightでツール作成

サンプルツール

簡単なサンプルツールをSilverlightで作成しましたので、こちらで最終形のイメージを掴んでください。(こちらがツールへのリンクになります。Internet exploreでご参照ください)

 Natural Cubic Spline Tool

上図のサンプルツールは、プロット上をクリックすると、該当するポイントのX-Yデータを読み取り、生成した各X-Yデータ間をNatural Cubic Splineで繋ぎます。

ただし、新しくプロットする点:\(x_{i+1}\)は、\(x_i < x_{i+1} \)である必要がありますので、ご注意ください。

上記のツールではNatural Cubic Spline補完をするため、生成した各Xのデータを8等分し、そのXにおけるYを補完で求めています。

ツールの実装

C#でも行列計算はMath.Net Numericsでも可能ですが、重くなるのでTDMAを入れ込みます。

当サイトに不具合、ご意見等ございましたらCEsolutionにお知らせください。

Page Top