OpenModelicaを使って電気回路のシミュレーションをしてみる
マルチドメインモデリング環境である、OpenModelicaを使用した電気回路のシミュレーションの例を通して、その使用感や言語的側面に触れたいと思います。
想定読者としてはModelicaに触れたことがない人で、バックグラウンドとして回路シミュレーションに多少なりとも興味がある人、という感じでしょうか。ちょっとした使い方に関するチュートリアルも兼ねているかもしれません。
Modelicaとは
オブジェクト指向のモデリング言語です。他分野に跨る複雑なシステムのモデリング(例えば、機械、電気、電子、油圧、熱、制御、等々)のシミュレーションを記述するのが得意で、微分方程式なども記述できます。
ざっくりいうと、常微分方程式で解けるくらいの数値解析のためのモデリングをするのに便利な言語です。また、適用分野ごとに使える標準ライブラリなども用意されており、手早くモデルを構築することができます。いろいろな分野の計算に使えるので、「マルチドメインモデリング言語」と言われることもあります。
Modelica言語の仕様は Modelica Association で策定されており、最近ではModelica 3.4が2017年5月にリリースされました。
Modelica言語を実装するシミュレータとしては下記のようなものがあります。(Wikipediaさんありがとう)
- CATIA Systems
- Dymola
- LMS AMESim
- JModelica
- MapleSim
- MathModelica
- OpenModelica
- SCICOS
- SimulationX
メジャーなシミュレーションツールである、MATLABがこの中に入っていないところは注目すべき点かもしれません。
OpenModelicaとは
Open Source Modelica Consortium(OSMC)が開発を進めている、オープンソースのModelica処理系です。OSMC は Modelica Association と連携して規格を策定しているようです。
OSMCはこの他にFMI(Functional Mockup Interface)規格というのを策定していて、規格の異なるモデル間の協調動作を規定しています。相互運用性を高めてMATLABを食っちまおう！という意気込みを感じるような気もしますが、どうなんでしょうか。
さて、ソフトウェアとしてのOpenModelicaに話を戻して、OpenModelicaの構成等について説明します。
OpenModelicaに含まれるもの
OpenModelicaは概ね、下記のツールから成り立ちます。
- OMEdit - モデルエディタです。この後、これを使用してモデリングを実施します。
- OMPlot - プロッティングツールです。シミュレーション結果を表示します。
- OMShell - Modelica言語のインタープリタを起動します。いわゆるREPLです。
- OMNotebook - OpenModelicaのドキュメンテーションのために作られたツールのようです。立ち上げるとOpenModelicaのクイックツアーが出てきます。
まぁ、OMEditが本体かな、というところです。
早速モデルを作ってみる
まずは立ち上げる
OMEditを起動します。
左側のペインにパッケージがたくさん並んでいるのが見えます。
モデルを作る
とりあえず "New Modelica Class" とかやって、モデルを作ります。そして左側のペインからドラッグ＆ドロップすると部品を落とし込めます。
部品のプロパティを設定できます。
出来上がったのがこちら。（3分クッキングかよ！）
昇圧チョッパー回路です。
感想としては、この手のツールのUIにしては使いやすい方です。
実行してみる。
とりあえず、実行。
とりあえず、何も考えずに一秒間、2000プロットくらいで試してみようと思います。
なにかエラーは出るが実行はされる。
適当に右側から項目を選ぶとプロットしてくれます。電圧を見てみます。
なんか、電圧跳ね上がってるし、落ち着いてないし、、、
電流も見てみます。
うーん、これはおそらくコンデンサのチャージ電流ですね。。。
(↑こいつのせいだ)
モデルをしげしげと眺めてみる
Modelica言語との関係
よく考えたら、Modelica言語全然出てきてない、って思います。(思うよね？)
ここまで表示してきたモデルは実はDiagram Viewというもので表示していました。これをText Viewに切り替えると。。。
Modelica言語による表現が出てきます。
なんとなく、HDL系の言語に雰囲気が似ていますかね。
これまで、ドロップしてきた部品のリストが上の方に出ています。実はペインの左側からドロップするという行為は、クラスのインスタンス化宣言をするということだったようです。また、部品のパラメータを入れるという行為はインスタンスの定数値を決定するということだったようです。
また、結線を行っていたものは下の方にconnect()という式によって表現されています。
なるほど。。。つまり、ここまで私はOMEditを使用してビジュアルプログラミングをしていたということになりそうです。
なお、annotation...となっている部分は表示上省略されていますが、これは図面の見え方を規定している模様です。
なんでこれで回路シミュレーションができるのか
さて、モデルを再度しげしげと見てみましょう。感想として「回路図やん！」というのが出てくると思います。（え？）
どう見ても回路図。。。
マルチドメイン言語じゃないの？回路だけ特別なの？って思ったりします。電流とか、電圧とか、それぞれ計算ルールが違うやん、それなのに、なんでつなぐだけで計算できるん？？
これ、ちょっとだけ説明しておこうと思います。
実は、電気部品はそれぞれでPinというコネクタをインスタンス化しています。Pinの定義は下記のようになっています。
VoltageとCurrentを持っているということです。あれ？Currentにはflowというキーワードが、、、
Modelicaの世界では通常、接続された物同士は等しい物として扱われます。PinコネクタにおけるVoltageでいうと、接続した点での電圧が等しいという制約を受けることになります。つまりこれはキルヒホッフの電圧則を表しているということになりますね。一方、flowを付与されて接続されたもの同士は和がゼロになるという制約を受けることになります。これは、、、キルヒホッフの電流則ですね。
こうやって接続された物同士に異なる制約を設けることによって電気回路の法則を表現しているわけです。
まとめ(個人の感想です)
OpenModelicaの特徴として下記が挙げられます。
- その辺の簡単な回路シミュレーター位のことはやってのけてしまいます。
- マルチドメインで異なるドメインをまたいだシミュレーションもできます。
- ユーザーインタフェースは比較的使いやすいと思います。(圧倒的とは言わん)
- 自分でModelica言語を書いて拡張することもできます。
- オープンソースなので問題があったらフィードバックすることも、自分で修正することも可能です。
- なにより、購入する手間もなく、ダウンロードしてすぐに使えます。
これは、使ってみない手はないのではないでしょうか。
おわび
スナップショットの撮り方が怠惰だったらしく、ウィンドウの影が残ってしまいました。ちょっと暗い四角が写り込んでいるように見えるのはそのせいです。すみません。
MathJaxの試し書き2
下記のお題で書く。積分や分数を含む式。
Q.169 ☆5
— どちゃ楽数学bot (@so_easy_math) 2017年7月26日
次の極限値を計算せよ。
Find the limiting value.pic.twitter.com/T9mOJmpvJN
十分大きな については
において
であるから
より
感想
ここまで来たら何で書いても多分しんどい。数学やってる人はえらい。
MathJax記法の試し書き
はてながMathJaxに対応しているらしいので、試しに使用してみる。
題材は下記の問題の解答を書くこと、としてみる。
Q.127 ☆8 [高2鉄緑校内模試]
— どちゃ楽数学bot (@so_easy_math) 2017年7月25日
(a+b)(b+c)(c+a)=2^2016を満たす自然数の3つ組(a,b,c)は幾つあるか。
ここでは、 に対して下記の一般解を求めることから始める。
として一般性を損なわない。
右辺の素因数は2のべきなので、左辺の各項（？）は
と書くことができる。ここで の大小関係から
が言える。
明らかに いずれかに奇数があるとすると、他の2つも必ず奇数となる。逆に、いずれかに偶数があると、他の2つも必ず偶数となることが容易に言える。
上記から は常に偶数であり、この差を
と表すと、(a)の第一式から
と書ける。また、これを(a)の第二式、第三式に適用すると
となり、
となる。
(b)より、 となるためには
が必要。(a)とあわせると
となる。
ここで奇数解の場合、 となりうるのは
の時のみ。すなわち、
となり、解は
となる。
一方、偶数解においては に対する解
がある場合、
に対する解 が必ず存在することが容易に確認できる。この関係を偶数解である限り繰り返したどることで、
に対する偶数解
に対して、
に対する奇数解
を得ることができる。逆にいうと、偶数解はある
に対する奇数解
を用いて
に対する解
と表すことができる。
前記の議論から、すべての解は奇数解 を用いて
と表すことができる。
(c)における奇数解の定式化と(d)における一般化を踏まえれば を用いて
により におけるすべての解が得られることがわかる。なお、対応する
は
であり、すべての は異なる解を生成する。
原題はの解の数を求めるのであるから、(e)より
を満たす
の組み合わせの個数を数えれば良い。
これに該当するのは の計336個である。
感想
- 綺麗にかけるけど、MathJaxって表示遅いな、、、
- タグはすべての環境で有効なんだろうか。不明。
- 細かいことだがインラインの数式の前後にスペースを入れるのを忘れてはいけない。HTMLそのものが西洋言語向けの文化に基づいちゃってるので仕方ない部分はあるんだろうなとは思うけどめんどくさい。
- TeXをよく知らないからなんだろうけど、displaystyleを使うときに前後に改行を入れるかどうかで左寄せかセンター寄せかが変わる？そういうもんなのか？
sh/bashのパイプ/リダイレクションの使い方
sh/bashのパイプ/リダイレクションの使い方についてよく使いがちなものを理解できるようにまとめておく。
モチベーション：挙動がよくわからない
シェルスクリプトに書かれているリダイレクションが理解できない。
これまで、運用とかしてこなかったので人の書いたスクリプトを読むこともなく生きてきて、シェルはtcshでできる部分でいいや、という感じで生きてきた。でも環境はいつもデフォルトで使っちゃう派なので何かしら自動化しようとするとsh/bashを使わないと色々めんどいことに気づいた。しかたがないのでいまさらながらリダイレクションの記法について学ぶことにしてみた。
自分の理解度テスト的に下記を上げてみたけど、こんな感じ。
cmd > file # わかる cmd >& file # なんだっけ？ cmd1 | cmd2 # わかる cmd1 |& cmd2 # わかる気もする。でも自分で書くときはこれが出てこない。 cmd 2>&1 > /dev/null # えっ？
以下に大体の読み方の基本から順番に書いていく。
基本
sh/bashのリダイレクションは基本的に (ディスクリプタ)> (ファイル) の羅列の形式である。
たとえば、
cmd 1> file
だとディスクリプタ1の出力先をfileという名前のファイルに設定する、という意味になる。
基本的にディスクリプタの1は標準出力、2は標準エラー出力なので、下記は「標準出力の出力先をfile1に、標準エラー出力の出力先をfile2に設定する」という意味になる。
cmd 1> file1 2> file2
>の代わりに>>を使うとファイルに追記する意味になる。
cmd 1>> file1 cmd >> file1
パイプは標準出力の出力先を後続のコマンドの標準入力に設定することを意味する。
下記はcmd1の標準出力をcmd2の標準入力に設定する。
cmd1 | cmd2
パイプにはディスクリプタ番号は指定できない。常に標準出力を次のコマンドの標準入力につなぐ。
ディスクリプタの複製
>& 記号を使ったディスクリプタの複製について説明するが、これは複製？よくわからないので複製というのは「そういう名前」だと思って深く考えないほうが読みやすいかもしれない。
a>&b はディスクリプタaの出力先をディスクリプタbと同じにすることを意味する。
下記は標準エラー出力を標準出力と同じターミナルに出力する。
cmd 2>&1
なお、ディスクリプタaは省略できてデフォルトは2である。bは省略不可。
つまり、下記も同じ意味になる。
cmd >&1
複数の記述がある場合の理解の仕方（これ重要）
基本パターン
複数の記述がある場合、基本的には左から解釈される。
この過程をこれから述べるように表の形式で理解していくと出力結果を考えやすい。
基本的に何もしなければ標準出力等の出力先は下記のように設定されている。
| ディスクリプタ | 出力先(当初) |
| 1 | 標準出力（ターミナル） |
| 2 | 標準エラー出力（ターミナル） |
ファイルへのリダイレクションが追加されると
cmd 1> file1 2> file2
| ディスクリプタ | 出力先(当初) | 1> file1 適用後 | 2> file2 適用後 |
| 1 | 標準出力（ターミナル） | file1 | file1 |
| 2 | 標準エラー出力（ターミナル） | 標準エラー出力（ターミナル） | file2 |
標準出力はfile1へ、標準エラー出力はfile2へ出ることになる。
複製を含む場合
複製を含む場合、上記で説明したように複製を「出力先を同じにする」と理解するのがポイントになる。
下記のような複製がなされた場合、複製は出力先を同じにするので下記のようになる。
cmd 2>&1
| ディスクリプタ | 出力先(当初) | 2>&1適用後 |
| 1 | 標準出力（ターミナル） | 標準出力（ターミナル） |
| 2 | 標準エラー出力（ターミナル） | 標準出力（ターミナル） |
下記のように順序を逆にすると結果が変わる。
cmd 2>&1 1> file1
| ディスクリプタ | 出力先(当初) | 2>&1適用後 | 1> file1 適用後 |
| 1 | 標準出力（ターミナル） | 標準出力（ターミナル） | file1 |
| 2 | 標準エラー出力（ターミナル） | 標準出力（ターミナル） | 標準出力（ターミナル） |
逆パターンだと
cmd 1> file1 2>&1
| ディスクリプタ | 出力先(当初) | 1> file1 適用後 | 2>&1適用後 |
| 1 | 標準出力（ターミナル） | file1 | file1 |
| 2 | 標準エラー出力（ターミナル） | 標準エラー出力（ターミナル） | file1 |
ディスクリプタの数を増やしてみる。新しいディスクリプタ3にはデフォルトでは何も設定されていない。
cmd 2>&1 3>&2
| ディスクリプタ | 出力先(当初) | 2>&1適用後 | 3>&2適用後 |
| 1 | 標準出力（ターミナル） | 標準出力（ターミナル） | 標準出力（ターミナル） |
| 2 | 標準エラー出力（ターミナル） | 標準出力（ターミナル） | 標準出力（ターミナル） |
| 3 | 未設定 | 未設定 | 標準出力（ターミナル） |
逆にしてみると
cmd 3>&2 2>&1
| ディスクリプタ | 出力先(当初) | 3>&2適用後 | 2>&1適用後 |
| 1 | 標準出力（ターミナル） | 標準出力（ターミナル） | 標準出力（ターミナル） |
| 2 | 標準エラー出力（ターミナル） | 標準エラー出力（ターミナル） | 標準出力（ターミナル） |
| 3 | 未設定 | 標準エラー出力（ターミナル） | 標準エラー出力（ターミナル） |
パイプを含む場合
パイプを含む場合、標準出力の当初値が変わる。
cmd1 | cmd2
| ディスクリプタ | 出力先(当初) |
| 1 | cmd2の標準入力 |
| 2 | 標準エラー出力（ターミナル） |
複製と併用するとこんな感じになる。
cmd1 2>&1 | cmd2
| ディスクリプタ | 出力先(当初) | 2>&1適用後 |
| 1 | cmd2の標準入力 | cmd2の標準入力 |
| 2 | 標準エラー出力（ターミナル） | cmd2の標準入力 |
この例は比較的記法の似た cmd 2>&1 1> file のパターンよりは cmd 1> file 2>&1 のパターンとの類似性を強く感じるものになっている。
よく「左から読め」と言われるのだけれども、パイプだけは右側に書かれているのに、一番左に書いたような挙動をするので非常にわかりにくい感じがする。
複製と他の構文をまとめた記法
複製には他の構文とまとめる記法がある。
ファイル出力とまとめて書くことができる。
cmd 2>& file
これは下記と同じ。
cmd 1> file 2>&1
パイプとまとめて書く記法がある。
cmd1 |& cmd2
これは下記と同じ。
cmd1 2>&1 | cmd2
復習
出だしにあったコマンド群
cmd > file # わかる cmd >& file # なんだっけ？ cmd1 | cmd2 # わかる cmd1 |& cmd2 # わかる気もする。書けない。 cmd 2>&1 > /dev/null # えっ？
は下記であることがわかるようになります。
cmd > file # 標準出力をfileに設定 cmd >& file # 標準出力と標準エラー出力をfileに設定 cmd1 | cmd2 # 標準出力をcmd2の標準入力に cmd1 |& cmd2 # 標準出力と標準エラー出力をcmd2の標準入力に cmd 2>&1 > /dev/null # 標準出力を/dev/nullに、標準エラー出力をターミナルに標準出力として出力
その他
2年ぶりに書いた記事がこれかよ、と悲しんでいる。
つらぽよいsecure boot
つらぽよAdvent Calendarにネタを提供したわけでも何でもないんだけど
自宅のデスクトップPCが壊れた。具体的に言うと、HDDがクラッシュした。
もうちょっと正確に言うと、動作真っ最中に子供がぶつかって倒してしまった。
昔のHDDと比べると、あのカツーン、って音がかなりささやかだけど、ありがちなあの音を立てて起動しない現象になった。
PC自体はLenovoのやっすいやつ(H520S)である。キーボードとか、スピーカーとか、こまいパーツコミコミで４万円とかしなかった。Core i5でWindows7が載ってた。
でも、我が家では結構ハイスペックマシンである。
どうしようか悩んだ。なぜか換装用のHDDはあった。SATA-II出だしの頃のやつで遅いけど、容量300Gとかあるので、デスクトップ機としては問題ない。Lenovoにはリカバリディスクを送ってくれるサービスがあるらしい。このHDDとリカバリディスクでハッピーになれそう。。。なわけなかった。7000円払う必要があるし、なによりLenovoの窓口自体が年末年始で休業中だ。
あー、じゃー、Linuxとかでいいじゃん。
Ubuntuのインストールディスクを手持ちのネットブックで焼いて、インストールしたら、無事インストールできた。さすが！余裕じゃないですか。いつもVMPlayerでやってるのと同じじゃん！
と思ってインストール後にPCを再起動したら、
Error1962:No operating system found.
ありー、なんかやり方まずかったかな。。。。
起動ディスクを別途作って起動しようとしても、全然起動しない。
ググってみると
Amazon.com: Customer Reviews: Lenovo H520s Desktop
Linux起動がうまく行かない、というような阿鼻叫喚の嵐。
集めた情報をだいたいかいつまんでまとめてみると
- 一部のLenovoのUEFI系マシンはsecure bootという仕組みがあり、Windows8が立ち上がるにはこれが必要である。
- 上記の仕組みを持つハードウェア環境でUbuntuなどの非正規OSはsecure bootにより起動を阻害される。RHELやFedoraの一部は正規の認証を受けているものがあるが、それ以外はsecure bootをOFFにしないと起動しない。
なるほど、secure bootをOFFにする必要があるのね、、、どうりで立ち上がらなかったわけだ。
と思い、BIOSのメニューを探してみたが、
secure bootに関するメニューすらない。。。。
マジか、、、同じ機種名の中でもOFFにできないモデルっぽい。
それでも、secureのまま、dual bootに成功している例はあるらしい。でも、逆に言うとdual bootの例しかない。どうも、Windowsのローダーを踏み台にして立ち上がるのは可能だが、それ以外はダメ、という感じらしい。(多分。たぶんね。)
つまり、WindowsのEFIローダーを書いておいて、起動シーケンスを横から頂戴すればいいのね！？
…いや、それやったら普通にリカバリディスク注文するし！
まー、どうしようかさんざん悩んだ挙句、Ubuntuのインストールディスクが起動していたのを思い出した。外部ブートはできるけど、光学ドライブのISOイメージだけってことっぽいんだよな。。。USBからも立ち上がんなかったし。
で、試しに、UbuntuのインストールディスクのgrubからHDD上のgrub.cfgを読み込んでみると、嘘のようにすんなりインストール後のHDDからシステムが起動した。
そういうわけで、当面、光学ドライブを踏み台にして起動する方向になりそうです。
なんだか変な感じではあります。
現時点ではUbuntuのインストールディスクを使っていますが、めんどいしアホいので、そのうちもっと小さなisoイメージ作ろうかなと思っています。
他にいい方法ないのかな。
正直、secure bootなるものにこういう風にいじめられるとは、なにか数奇な関係があるのかもしれません。
Evalモナドを使って無限リストの並列評価
http://partake.in/events/9f226986-2812-441d-98b7-1f5cca9be432 にオンライン参加してきたので、Evalモナドを使った並列評価を行なってみた。
コード
ナイーブな素数探索アルゴリズムです。指定された数以上の一定個数の素数を抽出します。
import Control.Parallel.Strategies -- rseq first n elements in parallel, then next n.. parListWithN :: Int -> [a] -> Eval () parListWithN n l = do parSeqListPar $ take n l case drop n l of [] -> return () _ -> parListWithN n $ drop n l -- rseq in parallel for finite length list parSeqListPar :: [a] -> Eval () parSeqListPar [] = return () parSeqListPar (hd:tl) = do rpar hd parSeqListPar tl rseq hd return () -- filter f l with parallel precalcuration for l filterEval :: (a -> Bool) -> [a] -> [a] filterEval f l = runEval $ do rpar $ runEval $ parListWithN 100 l return $ filter f l -- prime testing isPrime :: Integral b => b -> Bool isPrime n = not $ or $ map ((==0).(n `mod`)) [2..(n-1)] -- find 100 successive primes after 435300 main = print $ take 100 $ filterEvalN isPrime [435300..]
このコードの挙動のポイントはいくつかあります。
filterEval関数：並列処理のトリガーを引くコードを処理キューに積んで、全体評価を始めます。parListWithN関数：最初のn個の処理を並列処理にかけ、完了したら次のn個、というようにEvalモナドの挙動を指示する関数です。
並列処理がしたければ、rparを使ってタスクを積んでいくだけの簡単なお仕事です。
でも、次のような問題点があります。
- 全体評価の関数と並列処理が完全に分離していて、コントロールされていない。スケジューリング状況によっては、並列処理だけが進行して、全体評価側にフィードバックされずに進む可能性がある。
Evalモナドはこの辺りのコントロールが難しい感触があります。この後の並列計算ではこの辺りが解消されるか、気になるところです。
実行方法
使用時はthreadedオプションを付けてコンパイルします。
ghc -threaded FindPrime.hs
実行時は並列コアを指定します。
./FinePrime +RTS -N -RTS
threadscopeの実行結果とか見たかったけど、メモってる時点で時間切れなのでまた手が開けばやります。。。
モンゴメリ乗算の例
下記を参考に
モンゴメリ乗算 - Wikipediaを参考に、サンプルとして、モンゴメリ乗算の計算をしてみる。
`11*18(mod 169)`を計算する。
単純に、
11*18 = 198 = 1*169+29
で答えは`29`です。
これをあえてモンゴメリ乗算でやってみる。
モンゴメリ乗算は剰余算を減算回路の繰り返しなしに実現するテクニックです。剰余算、割り算が出てきますが、実際には2のべきしか使わないので、ビットマスクやシフト演算で高速にできるのが特徴。
`R`をひとつ決めると、サブの演算子を下記のようにとって
MR(x) = (x + (x*N' mod R)*N)/R Ninv : N' * N = -1 (mod R) となる数 R2 = R*R mod N
a * b = MR(MR(MR(a*R2)*MR(b*R2)))
とかけます。`R2`求めるときに`mod N` やってるやん！
R:256
に設定すると、下記のように関連する数値が導出されます。
Ninv:103 R2:133
`A=MR(a*R2)` `B=MR(b*R2)`とすると
A=(a R2+(a R2 Ninv mod R)N)/R
=(11*133+(11*133*103 mod 256)*169)/256
=112
B=(18*133+(18*133*103 mod 256)*169)/256
=45
AB=112*45=5040
MR(AB)=((AB)+(AB Ninv mod R)N)/N
= 157
MR(MR(AB)=(157+((157*103)`mod`256)*169)`div`256
=29あー、戻ってきた！
モンゴメリ乗算処理のテストケースを作ってみたのだけどそのままブログに書いてみましたの巻。