データ処理

回帰計算

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

NRTLのbinary parameterのRegression例

 

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

 

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

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

以下の流れでソルバーを設定します。

 SolverSet  

最終的なフィッティング結果です。

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

scipyでの回帰計算方法については、こちらをご参照ください。

 

初期値として、a12=1.0、a21=1.0、α=0.3を与え、a12~b21を推算します。(c12~e21とβ12は0で固定します。)

 

最終的なフィッティング結果です。

 parestPython  afterPython

回帰計算時の注意点

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

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

Data Reconciliation

最小二乗法と重みづけ

プラントの生データは、マテバラが取れていないことがままあります。

入出の流量計の誤差、分析精度の問題がありますので、シミュレーションモデル構築前にマテバラの取れたデータを用意する必要があります。

データ補完

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を使う場合は、逆行列が簡単に求まりますので、もっと簡単です。

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

PythonのInterpolate関数を用いる場合
 

こちらにまとめています。

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