はじめに
MATLAB,Python,Scilab,Julia比較するシリーズ。
といっても基本はベクトル行列ベースの演算に留め、
各環境独自の機能はあまり使わない方針。
あくまでベクトル＆行列でどこまでできて、
ベクトル＆行列に留めているが故に各ツール/各言語でやっても同じことができそうだ。
ってことを確認することが目的。
ちなみに、問答方式じゃない方も記事もある。
問答方式に合わせて内容を強化していく予定。
書籍とか
その他の章
概要偏
- ベクトル、行列良く分からん。
- 単なる表現ルールではあるが、そのルールがわからん。
- 学生時代は大量の手計算をさせられた記憶がフラッシュバック。
- 今の時代はツール使えば計算は一撃。
- MATLAB,Python,Scilab,Julia並行でいろいろやるプランを考える。
- まずはMATLABの概要説明。
- 自動車業界だとおなじみのツール。
- 価格はお高め。
- homeライセンスというのもあり、こちらはかなりお安い。
- ただし、coder系は含まれていない。
- homeライセンスというのもあり、こちらはかなりお安い。
- Pythonの概要。
- PythonというかNumpyがベクトル、行列に強い。
- for文で回すよりベクトル、行列として演算させた方が圧倒的に速くなることは多い。
- Python自体は無料と思って良いが、パッケージ別にライセンスが規定されてるので要注意。
- ScilabはMTALAB/Simulinkと比較されることが多い。
- 理由は、xcosというグラフィカルツールにてブロック線図を書けるツールの存在。
- 互換性は全く無いが、Simulinkっぽいツール。
- MATLAB/Simulinkを自宅で使えない場合に使ってる人は多かったのでは？
- ぶっちゃけJuliaについては良く知らない。
Wikipadiaに活躍してもらって引用しまくり。 - 最大の特徴はJITコンパイル。
- これにより、複数回の実行に関しては類似言語/ツールからみたら最速になり得る。
- というわけでJuliaに関しては逐次調べながらの試しながらで進める。
行列
- そもそも行列の存在意義に疑問が・・・。
- 分かり易い存在意義としては連立方程式がわかりやすい。
- 行列を使うとルールベースで解くことが可能。
- ここらへんについては次回説明予定。
- とりあえず連立方程式を普通に解いてみた。
- しかし、そのプロセスをプログラム化するのは超難解。
- つまり、行列を使うとこの超難解な状態から脱することができる。
- 今回は、連立方程式を行列で表現するとどうなるか。ってところまで。
- 行列で連立方程式を解いてみた。
- 両辺を行列で割る・・・のがだが、行列は除算が無く、逆行列を掛けるで除算を実現する。
- 逆行列は掃き出し法と呼ばれる方法で求めるが、ここでは公式を使用。
- 結果として、答えが求まった。
- 特徴としてはルールが明確なためプログラム化し易いという点。
- 逆行列は掃き出し法にて求めることができる。
- 実際に掃き出し法を実施。
- 前回使用した逆行列が求まった。
- MATLAB、Python(Numpy)、Scilab、Juliaでは逆行列を求める機能があるので、直に計算することはない・・・想定。
- 連立方程式を解くということは複数の関数の交点を求めるということ。
- 行列はそれを一撃で解ける。
- ためしにMATLABで算出したら一撃。
- 移動体の予測線を関数と見なすと、交点を求める重要性がわかりやすいかも？
基本的な使い方
- 各ツール、言語の基本的な使い方として以下をやってみる。
- 単純なスカラー計算。
- ベクトルの定義。
- 等差数列の作成。
- 行列の定義。
- まずは手馴れたMATLABで実施。
- Python(Numpy)とScilabの基本的な使い方。
- Python(Numpy)は以前から使っている物なので手馴れたもん。
- ScilabはMATLABと同一の記載方法でいける。
- ただし、コメントアウトが「%」じゃなくて「//」
- ここも一緒だと楽だったが・・・。
- ただし、コメントアウトが「%」じゃなくて「//」
- Juliaの基本的な使い方・・・の前にいろいろクセが違うのでそれの調査。
- start:step:endの形式(区間演算子)で等差数列を表現できるが、この状態ではメモリ上に実態を持っていない。
- よって、読み出しはできるが、書き込みはできない。
- 区間演算子に実態を持たせるにはVectorに渡すことで解決。
- Juliaの基本的な使い方。
- Juliaは列ベクトルがデフォルト。
- MATLAB、Scilabは行ベクトルがデフォルトであるため、扱いに気を付ける必要がある。
- 列ベクトルがデフォルトになっている理由としては、数式との一致性を考慮した結果と推測される。
- 基本的な使い方の続きとしてスライシングについて。
- 特定の要素、特定範囲を抽出可能。
- 区間演算子start:step:endを元に範囲抽出するが、step=1なことがほとんどなので、stepを省略したstart:endの書き方になることが多い。
- Python(Numpy)でスライシングを実施。
- 0オリジンのためMATLABと設定する数値が異なる。
- 加えて、区間演算子の終端は範囲に指定範囲には含まれない点に注意。
- Scilabでスライシング。
- MATLABと同一。
- Juliaでスライシングを実施。
- 基本的にはMATLABと似た感じ。
- ただし、配列添え字用のカッコが違う。
- あと、スライシングの結果、ベクトルとなった場合は列ベクトルになる。
- 行列としてスライシングした場合は、元の行と列の関係は維持される。
行列演算
- 代表的な行列演算を列挙。
- 基本的な四則演算に加えて、アダマール積、べき乗、転置。
- まずは加算、減算。
- 各要素単位で加算、減算すればOK。
- 当然、「次元を一致させる」必要がある。
- 行列の乗算(内積)について説明。
- 上記はなぜそのような演算になるか不明(太郎くん談)。
- これを理解するには線形代数の基礎部分を理解する必要がある。
- 線形代数すべてを説明するとなると大変だが、基礎部分を可能な限り簡単に説明予定。
- 今回はアダマール積について。
- 演算子がいろいろあり、アダマール積かどうかは文脈で読み解くしかない。
- しかし、特殊な状況でしか登場しないので、そういうものがあるという程度で留めておいてもよいかも。
- 画像処理の畳み込みで出てくることは多い。
- 行列の除算について。
- 行列は原則的に除算は存在しないが、「逆行列を掛ける」がそれに該当する。
- さらに行列の積は結合法則はあれど、交換法則はない。
- 上記に伴い、左除算、右除算と言う概念が出てくる。
- 逆行列の位置が変わる。
- 数式上ではあまり出て来ないが、各ツール、言語がサポートしていることが多い。
- 行列の転置について説明。
- 転置自体は、行列の行と列を入れ替えるだけの話。
- 具体的な利用シーンというのは特になく、計算都合で使うことがほとんど。
- 良く使う処理なので、名前が付いていた方が利便性が良いという考え方が妥当そう。
- やっぱり線形代数の基礎はやっておく。
- 割とすぐに詰む可能性があるから。
- 説明手順をとりあえず決めた。
- 行列の内積の公式の再確認。
- 方程式と内積。
- 連立方程式と行列。
- 行列によるベクトル変換。
- 行列によるベクトル群変換。
- 行列の内積の公式の再確認。
- 行列の内積の公式の再確認。
一旦忘れてOK。
- 内積の定義を確認。
- 内積は単なる計算方法であり、内積そのものにに意味はない。
- ただし、特性のようなものはある。
- 内積の分かり易い特性としては相関性。
- 類似度とも言われ、特に内積を利用したものをcos類似度と呼ばれる。
- 基本的な計算であるが故に畳み込み積分、類似成分抽出、方程式などに利用される。
- 余弦定理を何とか証明。
- 垂直線を使って2つの直角三角形を作ることで各辺を三角関数を使用した表現が可能。
- 三角比の基本公式を加えると、余弦定理が求まる。
- 基本公式は三平方の定理と半径1の円起動の点と原点をを元に作った直角三角形から求まる。
- 内積の定義と余弦定理から成分表記の内積を求めた。
- ベクトルとしての内積と、成分表記としての内積が等しいことを証明。
- 上記を利用して、内積が方程式と強い関係性があることを示すの次回。
- ベクトル内積の公式を再掲。
- ベクトルの内積で方程式を表現できる。
- n次方程式、多変数方程式でも考え方は一緒。
- 二元一次方程式を書き出す。
- 上記をベクトルの内積で表現し直す。
- さらに上記を行列で表現し直す。
- まずはこれが最もシンプルな行列の性質を示している。
- 連立方程式を行列演算にしたものを再掲。
- 上記の構成は[出力ベクトル]=[変換行列][入力ベクトル]となる。
- これの入力、出力を列ベクトル2セットにすると2×2の行列になる。
- すべて2×2行列になるが、数式上の位置によって、入力、変換、出力と意味が異なる。
- このルールをすっ飛ばしてる行列嫌いになるかも？
- 前回の行列演算の式を再掲。
- 前回は入力、出力の列ベクトルを2セットだったが、これをnセットにするとどうなるか。
- 入力、出力が2×2行列からnx2行列へ。
- これの利点は入力から出力の変換が一括で表現できること。
- これを知ってるだけでコミュニケーションスキルが上がる？
- これの利点は入力から出力の変換が一括で表現できること。
- 入力、出力が2×2行列からnx2行列へ。
- 具体的に各ツール、各言語での行列演算についての話に突入。
- 四則演算に加えて、アダマール積、べき乗、転置、反転。
- まずは手馴れたMATLABから。
- Python(Numpy)での基本的な行列演算を確認。
- 内積は「@」
- 「*」だとアダマール積になるので注意
- それ以外にも0オリジンだったり、終端指定が-1だったりとクセが違う。
- 右除算もできるように見せかけて、実は逆行列とのアダマール積なので目的としたものとは異なる。
- Scilabで基本的な行列演算実施。
- 基本的にはMATLABと同一。
- ただし、要素終端がendではなく「$」である点に注意が必要。
- Juliaで基本的な行列演算を実施。
- 大体MATLABと一緒だが、以下の違いがある。
- 配列添え字のカッコが丸カッコじゃなくて角カッコ。
- flipdimは使えなくて、代わりにreverseという関数を使用する。
- 以前はflipdimは存在していたようだが、現在では無くなってる。
ユーザ関数
- ベクトル、行列から離れて、少しプログラミングより話にシフト。
- 各ツール、各言語でユーザ関数の作成方法を確認する。
- MATLABは関数名と同名のmスクリプトファイル名にする必要あり。
- Pythonの場合のユーザ関数作成方法。
- 対話モードで作成する場合とスクリプトファイル上で作成する場合がある。
- が、実際は対話モード時のルールが共通で適用されてるだけ。
- 他のスクリプトファイルで定義した場合はimportを使用。
- エイリアスで名称変更可能。
- 対話モードで作成する場合とスクリプトファイル上で作成する場合がある。
- Scilabの場合のユーザ関数作成方法について。
- MATLABと似ていると思いきや、全く異なる仕組みっぽい。
- スクリプトに記載したとしても、明示的にワークスペースに関数を展開する必要がある。
- 仕組みは異なるが、関数として展開してしまえば使い方は一緒と言える。
- Juliaの場合の場合のユーザ関数作成方法について実施。
- 基本的にはPythonに似ている。
- 2変数以上を戻す場合は、明示的にreturn文を使用する必要がある。
- 他のファイルで関数を定義している場合はinclude文を使用する。
- C言語のincludeに似ている。
- ユーザ関数定義に引き続き、波形表示も人間向け機能。
- MATLABによる波形表示を確認。
- plotで表示。
- hold onで同一グラフに表示させる設定が可能。
- subplotでグラフ分割。
- ラインスタイル、色、マーカの指定が可能。
波形表示方法
- pythonで波形表示する場合はmatplotlibを使用する。
- matplotlibはMATLAB仕様に合わせこんでくれている。
- マーカに関しては、MATLABにはない指定子もある。
- Scilabの波形表示はMATLABと同一。
- 特殊なグラフ表示は乖離する可能性が高いが、そこまで複雑使い方はしない予定。
- JuliaはPlotsかPyPlotで波形表示。
- PyPlotはmatplotlibのラッパーらしく、使い勝手が他の環境と似ている。
- PyPlotの描画パラメータは個別に指定する必要あり。
状態空間モデル
- 状態空間モデルの話に突入予定。
- その前に微分の記法について疑問点浮上。
- 微分記法は以下がある。
- ライプニッツ記法。
- ラグランジュ記法。
- オイラー記法。
- ニュートン記法。
- 状態空間モデルではニュートン記法が一般的。
- 暗黙的に時間微分であることがわかるため。
- 状態空間モデルではニュートン記法が一般的。
- 状態空間モデルの数式について説明。
- 方程式。
- 状態方程式と出力方程式。
- 変数。
- 状態量、入力量、出力量
- パラメータ。
- 状態行列、入力行列、出力行列、直達行列。
- 方程式。
- 各呼び名が揺れるのは使用する領域が広く、観点が異なるため？
- 状態空間モデルの各要素は分かれど使い方はわからない。
- 使い方を見てもよくわからない。
- よって、超シンプルな微分方程式を対象に状態空間モデルを作ってみる。
- ニュートンの運動方程式を対象とする。
- まずは状態量を定義。
- 速度、距離を状態量とした。
- 運動方程式を紐解く。
- 距離、速度、加速度の関係性が微分を吸収する。
- 状態量の内訳である速度、距離の方程式が求まったところ。
- 状態方程式、出力方程式を組み上げた。
- 状態方程式は前回の運動方程式から導出した微分方程式を元に作成。
- 出力方程式は参照したい状態量に合わせて出力行列Bを調整するのみ。
- 状態空間モデルを確認するにはシミュレーションしてみるしかない。
- まじめにシミュレーションしようと思うとベクトル、行列に対する微分を解決する必要がある。
- (これもやる予定だが後で)
- 各ツール、各言語で状態空間モデルが扱えそうなので、それらで動かしてみる。
- ただし、MATLABに関してはSimulinkの状態空間モデルブロックで実施予定。
- MATLAB/Simulinkで状態空間モデルのシミュレーション。
- 必要ブロックはState-Space、Step、Scope、Mux。
- 状態空間モデルの各行列の設定はState-Spaceの詳細設定で可能。
- シミュレーションは摩擦等を無視しているので宇宙空間での挙動と思って。
- Pythonで状態空間モデルを扱うには、controlライブラリのmatlabモジュールが必要。
- 仕様的にはMATLABのControl System Toolboxを踏襲している。
- ss関数に各行列を渡し、システムオブジェクトを取得。
- lsimに入力のstep信号をシステムオブジェクトを渡してシミュレーション。
- Scilabで状態空間モデルのシミュレーションをするにはsyslinとcsim関数を使用する。
- MATLABに寄せてるかと思いきや、この分野はかなり異なる仕様になっている。
- 想定通りのシミュレーション結果を得られた。
- Juliaで状態空間モデルをシミュレーション。
- Pythonと同じくMATLAB Control System Toolboxの仕様を踏襲したControlSystemsパッケージを使用。
- 他のツール、言語と同じ結果が得られた。
- 状態空間モデルを掘り下げる。
- ベクトル、行列であることの利点があるはず。
- その前に状態空間モデルを使用しない場合はどうなるかをやってみる。
- 漸化式、ブロック図、C言語化をやってみる。(漸化式は今回やった)
- 運動方程式をブロック図で書き起こした。
- エディタの利便性からSimulinkを利用。
- 上記を離散化。
- Simulink利用の場合、この段階からオートコード生成。
- 離散化ブロックを元に漸化式。
- これからCコード化をする。(次回Cコードを見せる予定)
- ブロック図、漸化式を再掲。
- 上記を元にCコード化。
- 上記のシミュレーション結果を見せた
- 状態空間モデルのまま微分解決可能なはず。
- まずは状態方程式の微分解決を実施。
- 両辺を積分して、その後にオイラー法で微分解決。
- とりえずはオイラー法でも精度が十分なことは多い。
- 両辺を積分して、その後にオイラー法で微分解決。
- 出力方程式の微分解決を実施。
- とはいっても、出力方程式側には微分方程式は居ないのでxを代入しただけ。
- これで、状態空間モデルのままで演算できる状態になったと言える。
- しかし、これでもプログラム化のイメージは湧きにくい。
- よって、次回から各ツール、各言語で書くとどうなるかを確認。
- MATLABでベクトル、行列演算による状態空間モデルの演算実施。
- 導出した数式のまんまでコードが組める。
- このルールに即していれば、さまざまな振る舞いを規定できる。
- シミュレーション結果も想定通り。
- Python(Numpy)でベクトル、行列演算による状態空間モデルの演算実施。
- 流れとしてはMATLABと同一。
- 内積の演算子が「@」な点に注意。
- シミュレーション結果も想定通り。
- Scilabでベクトル、行列演算による状態空間モデルの演算実施。
- MATLABと同一。
- グラフ表示の部分に差異があるだけ。
- シミュレーション結果も想定通り。
- Juliaでベクトル、行列演算による状態空間モデルの演算実施。
- MATLABとほぼ同一。
- 添え字、ドット演算子に違いあり。
- シミュレーション結果も想定通り。
- Cコードによるベクトル、行列演算による状態空間モデルの演算実施。
- MATLAB Coderで出力。
- シミュレーション結果も想定通り。
- コード自体は複雑になるが、多変量の微分方程式になった際に効果は大きくなる。
- 状態空間モデルでもうちょっと複雑なものを取り扱う。
- とりあえずDCモータを採用。
- 導出の流れは運動方程式の時と一緒。
- まずは状態量を特定。
- 状態量の特定は入力から期待する出力に至るために必要なパラメータをイメージすると良い。
- 各種必要な微分方程式を特定
- 角度→角速度、電流→角速度、を特定。
- 電圧→電流はオームの法則とキルヒホッフの第2法則の組み合わせで導出する。
- とりあえず、電気回路として見た場合と、逆起電力の部分を特定。
- 次回これらを合体させる予定。
- とりあえず、電気回路として見た場合と、逆起電力の部分を特定。
- キルヒホッフの第2法則に則って電圧導出式を合体さえた。
- 単純に加算。
- しかし、今回は逆起電力なので、結果的には引き算にはなる。
- 単純に加算。
- 電流から電圧を求める式に変形。
- 電流が1階微分されているが、元々電流の1回微分が欲しいのでちょうど良い。
- 状態量と各種微分方程式を再掲。
- 上記の情報から状態方程式を組み上げた。
- 表現がベクトル行列になっただけで、導出した微分方程式と一緒。
- 状態空間モデルも見慣れてしまえばそれほど不可思議なものではない。
- 出力方程式を導出。
- 基本は、出力させたい項を出力行列で指定するだけ。
- 実験段階では全状態を見たいことが多いので全部出力指定にすることが多い。
- 各種行列を列挙。
- シミュレーションしないと確からしさはわからないが、以前のシミュレーションコード使い回しでいけそう。
- これが状態空間モデルの良いところ。
- DCモータ状態空間モデルをMATLABでシミュレーション。
- 状態空間モデルを演算する関数自体はそのまま使い回し。
- シミュレーションとしては想定通りの結果。
- DCモータ状態空間モデルをPython(Numpy)でシミュレーション。
- 流れとしてはMATLABと一緒。
- 状態空間モデルの演算用関数が変化しない特徴も一緒。
- シミュレーションも同一であり、想定通り。
- DCモータ状態空間モデルをScilabでシミュレーション。
- 演算自体はMATLABと同一。
- 差はグラフ表示の微調整のところ。
- シミュレーションも同一であり、想定通り。
- DCモータ状態空間モデルをJuliaでシミュレーション。
- かなりMATLABと近似のコードになる。
- linspaceがrangeになってるくらい。
- シミュレーションも同一であり、想定通り。
状態空間モデルをPID制御
- DCモータ状態空間モデルを制御すべくPID制御器を追加予定。
- 上記に至るロードマップ提示。
- 「PID制御器の離散化」がやや数式まみれになるためちょっと鬼門。
- オイラー法による微分解法するだけなので比較的簡単にするつもり。
- これまでの状態空間モデルの構成を確認。
- オープンループ制御になっている。
- 今回の構成を確認。
- PID制御でクローズループにする。
- 角速度ωをフィードバックすることで角速度制御をすることが分かる。
- PID制御の基本式を確認。
- 上記の積分を外側に追いやった変形式を確認。
- 積分を外側に追いやる方がワインドアップ対策がし易いので、こちらが使用されることが多い。
- PID制御の数式をオイラー法で微分解決して離散化。
- 一見するとカオスには見えるが、似たような式が並んではいる。
- 上記がブロック線図にする際にいい感じに効いてくる。(たぶん)
- PID制御離散化式からブロック線図を作成。
- 同一の数式が埋まっている箇所は同一の信号線を参照できる。
- よって、比較的にシンプルになる。
- 同一の数式が埋まっている箇所は同一の信号線を参照できる。
- 正しさの証明はやはりシミュレーション。
- 構成図、ブロック線図、ソースコードの3点セットで見て行けばそれほど混乱はしない。(たぶん)
- PID制御器のブロック線図と全体構成を再掲。
- 上記をMATLABで実現。
- 接続に関してはコード上では分かり難いので全体構成図と見比べながら確認した方が良い。
- シミュレーション実施。
- u(t)の挙動と見るとPID制御っぽい挙動になっている。
- 各PIDゲインを調整すると挙動が変わるはず。
- u(t)の挙動と見るとPID制御っぽい挙動になっている。
- MATLABでやったDCモータ状態空間モデルをPID制御をPython(Numpy)で実施。
- Pythonの場合、構造体はclassで実現。
- 事前にclassを定義する必要はある。
- MATLABと同様の結果が得られた。
- ScilabでDCモータ状態空間モデルをPID制御シミュレーション実施。
- コード自体はMATLABと一緒。
- 構造体の生成の仕方も一緒だが、内部的には連想配列で実現されている。
- シミュレーション結果もOK。
- JuliaでDCモータ状態空間モデルをPID制御シミュレーション。
- Juliaも構造体を定義できるが、structだと更新不可になる。よってmutableを使って更新可能な構造体にする必要がある。
- シミュレーション結果は他のツール言語と同じ結果が得られた。
コメント