役立つTipsを紹介します。コンテンツは少ないですが、Scipy、Numpy、Matplotlibの基本的な使い方に関するTipsを順次追加しています。※2017.06.10追記:なにはさておきAnacondaとPycharmをインストールしてください。numpyやscipyのインストールで悩む必要がなくなります。個人的にPythonが有難いのは以下の点です。1行列計算や多次元配列の取扱いが簡単(numpy)2Scipyの豊富な計算ライブラリ(最適化関数やODEのSolver等)3グラフ描画が比較的容易で綺麗に描ける(Matplotlib)4複素数計算が可能(Bode線図やNyquist線図を作る際に便利)5外部ファイルの入出力が行いやすい6変数宣言をせずに使える(コードが短く書ける)Pythonでの数値計算の恩恵を受けるためにいくつか便利なライブラリを以下に紹介します。・SciPy、NumPy、matpoltlib:最適化計算や逆行列演算、グラフ描画等に便利なライブラリです。・pandas:データ処理に便利なライブラリです。・PySide:GUIの作成に便利なライブラリです。ちなみに私はインストールで何度もはじかれて苦戦しました。インストールについてはゆくゆく整理していきたいと思います。※2016.05.08追記:Scipy、Numpy、Matplotlibのインストールについては、こちらがよくまとまっています。※Numpyのインストールが必要です。\end{bmatrix}$$a=np.array([[1],[2],[3]])また上記の行列を転置した以下の行列はa=np.transpose(a)のように書けます。\end{bmatrix}$$行列同士の積は、例えば行列aと行列bの積であれば、np.dot(a,b)と書くだけです。\end{bmatrix}\)と、その転置行列の積をPythonで計算する場合、非常に短いコードで済みます。出力例このようにNumpyを使うと、転置行列の作成、行列の積が簡単にできます。・csvファイルの読み込みについてnumpyのloadtxtを使うと簡単に読み込むことができます。np.loadtxt("sampledata.csv",delimiter=",")numpyのloadtxtを使うと簡単に読み込むことができます。csvファイルを読み込んで、データの統計処理を行うのも簡単です。最大値print(np.amax(a))最小値print(np.amin(a))算術平均print(np.mean(a))中央値print(np.median(a))標準偏差print(np.std(a))分散print(np.var(a))Matplotlib基本コードこちらのサイトが完璧です。私がよく使う機能だけを以下まとめておきます。"blue"色の種類はこちらにまとめられています。基本的な色は以下テーブルのように、そのまま色の名称を記載するだけです。以下にdark、lightをつけることで色の深みを変更することも可能です。shadow=True)こちらとこちら、こちらが参考になります。Gridの表示こちらが参考になります。・2dplot出力について詳細はゆくゆくまとめていきたいと思いますが、matplotlibのplotはsubplotをまとめて一つの形式で出力できるため、非常に便利です。NRTLで水-EthanolのXY、PXY、活量係数をまとめてプロットしたいという場合、以下のようにまとめて表示させることが可能です。以下のように各subplotの配置位置を2dgridで記載し、各suubplotの詳細を記載していきます。・3dplot出力についてQiitaのブログに詳細記載しましたのでご参照ください。例として、以下のように理想気体と実在気体のPVT図を重ね書きする方法をまとめています。pressureのデータ(データはperryのhandbook参照)をテキストファイルに保存しておき、各係数を必要に応じて読みだすことを考えます。$$ln(P)=C1+C2/T+C3ln(T)+C4T^{C5}$$C1:73.649、C2:-7,258.2、C3:−7.3037、C4:4.1653E−06、C5:2ただし、圧力の単位はPa、適用温度範囲は、273.16~647.096[K]となっています。上述のC1~C5パラメーターを各行毎に区切ってテキストに入力しておけば以下のような形で読み込むことができます。1f=open('VPwat.txt')で同階層にあるVPwat.txtを開き、2lines2=f.readlines()でテキストファイルを1行ずつlines2に格納していきます。※readlines()以外にも、read()、readline()があります。使い分けは以下のリンクが分かりやすいので、ご参照ください。read()1f.close()でテキストファイルを閉じます。4他の蒸気圧推算法でも対応できるように、VPの係数はMax7個まで格納できるようにしています。この際、Cに格納するデータが空だとエラーが出るため、空の部分には0を代入しています。5C1~C7までデータを格納していきます。この際、リストからそのままデータを引っ張ってくると改行コードが邪魔であるのと、数値として認識されなかったため、.strip()で改行コードを除き、float()でそれぞれを数値化しています。このように、係数は分離して別置きしておけば、推算式が変わっても呼び出し元の関数を使い分けるだけで対応することができます。・物性基本情報を記載したExcelファイルの読み込みについてExcelファイルの読み取りまずxlrdを読み込みます。xlrd例えば、物性データを「PhysicalData.xlsx」として保存し、ここから必要な情報を読み取る場合、以下のように設定します。ExcelデータPythonコード簡単に各行のコードをご紹介します。・interpolate(補間)こちらのリンクに使い方がまとめられています。Akima1DInterpolator(x,y)こちらのリンクで紹介されている関数:\(sin(20x)+exp(5x/2)\)について各補間法を試してみます。・optimize(leastsq関数)こちらのリンクを参考にさせていただきました。1次遅れプロセスのシステム同定1次遅れプロセスのシステム同定をleastsq関数で求めます。同定する1次遅れプロセスを、以下の伝達関数と仮定し、ステップ信号を与えた場合の時間応答データから、各パラメーターK、τ、Cを求めていきます。パラメーターの初期値)の形で用います。適当なパラメーターについて、伝達関数のステップ時間応答のデータをツールで出力します。以下、K=1、τ=5、C=1とした場合の出力結果です。出力結果を一旦エクセルにコピペし、乱数を与えて数値を少しズラした以下のデータをサンプルデータとして与えます。(csv形式で保存し、配列にデータを読み込みます)初期値を、K=1、τ=1、C=1として、各パラメーターをleastsq関数で求めたプロットが以下になります。データのバラつきをあまり強くしなかったので、サンプル点と殆ど重なって見えていますが、推算パラメーターは問題無さそうです。leastsq関数で求めた各パラメーターは、以下の通り。983.033708335プロットからも明らかですが、本来の解:K=1、C=1、τ=5とほぼ同等の値が得られており、時定数における出力値が1-1/eとなっていることが確認できますNRTLのバイナリパラメーター回帰もう少し複雑な関数の最適化事例としてNRTLのパラメーター推算をこちらに整理しています。・odeint,ode(常微分方程式のsolver)こちらのリンクを参考にさせていただきました。上記のリンクによると、odeinitを使うと、硬い問題かどうかによって、採用手法を自動で切り替えてくれるとのこと。odeを使うと、解法を自由に選択できるようなので、まずはodeを試してみます。主要な解法を以下に記載します。method='adams')lsodaの呼出し:ode(function).set_integrator('lsoda')ode45の呼出し:ode(function).set_integrator('dopri5')(functionの部分に関数定義します。Pol方程式で各手法を比較してみます。pol方程式は以下の通り。\left(1-x^2\right)\frac{dx}{dt}+x=0$$Adamsは今回の条件では解くことがでなかったため、その他の計算手法の計算速度を比較してみると以下のようになります。BDF:19.71[sec]、lsoda:0.14[sec]、ODE45:13.73続いてodeintを試してみます。t)の形で定義します。(X0が初期値、tが時間定義になります)計算時間わずか0.0164[sec]。基本的にはまずodeintを使うという戦略は妥当そうです。上記VDPの解をmatplotlibで作図すると以下のようになります(BDF)。・signal(Bode線図の作成、時間応答)こちらのリンクを参考にさせていただきました。phase)Plot結果は以下の通り、正しい結果が得られています。時間応答(Impulse、Step関数)time,y=signal.impulse(伝達関数)time,y=signal.step(伝達関数)Plot結果は以下の通り、正しい結果が得られています。Technic.こちらのサイトが詳しいです。Regressionこちらのサイトを参考に、無駄時間1次遅れプロセスの表現が可能か試してみました。Excelで1次遅れ&無駄時間プロセスの時間データを与えて、SVRで確認してみます。カーネル関数はRBFとPolynominalを試し、上記リンク同様、以下のパラメーターを採用しています。2)雰囲気は出ていますが外挿範囲も含め、さすがに無理がありますね。・もう少しトレースしやすそうな例としてScipyの活用でご紹介した1次遅れプロセスの表現をRBFでも行ってみます。例1)1次遅れプロセスの表現こちらのサイトを参考にさせていただきました。・1次遅れプロセスの時間データを学習の入力値とし、1次遅れ時間、プロセスゲインを出力値として学習させます。どのくらいのデータが学習に必要か分からないので、1次遅れ時間を変更したケースを10ケース、ゲインを変更したケースを10ケース、それぞれについての時間データを20個用意し、学習ケースとして用意しました。データの生成は以下のようにエクセルで行いました。(for文で回しながら順番にcsv出力)各csvデータを入力として書き込みます。例2)無駄時間1次遅れプロセスの表現おまけ・絶対に勝てない石取りゲーム息子用に遊びで作ったコードです。意味は無いですが、対話型コードの参考として見てみてください。While,If文の用法も確認でき、意外とこうした遊び感覚で作るプログラムの方が勉強になったりします。Shellでの入力要求はinput()の形で書きますが、そのままでは文字列として認識されてしまいます。このため、数値として使用するためには、int(input())として一旦数値に変換する必要があります。逆にprintで出力するには、文字列に変換する必要があるため、str()を使用します。当サイトに不具合、ご意見等ございましたらCEsolutionにお知らせください。