Modelicaというと,LEGOのように色々な機械要素のコンポーネントを組み立ててモデルを作っていく「ソフトウェア」だと思われることも多いのですが,ModelicaはFORTRANやC++やPythonと同じように「プログラミング言語」の一種です.ただし一般的なプログラミング言語とは異なる側面も持っており,物理現象を記述するのにより適した仕様になっているため「モデリング言語」とも呼ばれています.現在は商用のものも含めて多くの1DモデリングツールがこのModelica言語をベースに作られています.1Dモデリングに取り組もうとする方はまず,Modelicaとはソフトウェアではなく言語であるという事を理解する必要があります.
スクリプトベースとは?
Modelica言語をベースにModelica協会によって開発されたオープンソースの1DツールがOpenModelicaです.OpenModelicaだけでなく,ほとんどの1Dモデリングツールは下記のようなグラフィカル・インターフェースを搭載しています.企業においてシステマティックなシミュレーション環境を構築していくためには,生産性とメンテナンス性に優れたグラフィカルベースによる1Dモデリングはやはり必須と言えるでしょう.
しかしながら(オープンソースのOpenModelicaであっても商用ソフトであっても)基本的なスクリプトでの理解を飛ばしていきなりグラフィカルベースから入ってしまうと(例題のような単純なモデルであればいいのですが)実務で通用するような応用的なモデルの構築は逆に遠回りになります.下記は上記のモデルをスクリプト表示させたものです.
という訳で,まずは基本的なスクリプト・ベースによる1Dモデリングを理解するところから始めましょう.ここではシンプルな機械系の1Dモデル(マスバネダンパーモデル)を使って説明していきます.
作成する工学モデル
下記は,機械のダイナミクスを表現する基本的な支配方程式です.実際の機器を扱う場合はここにさらに要素が追加される事も多いですが(例えばクーロン摩擦項とかストライベック効果項とか)基本的な動力学の構成要素である慣性力,減衰力,復元力(剛性)は,下記の数学モデル(運動方程式)で表現する事ができます.
$$\large{
\begin{array}{r}~~~\underline{M~a}~~~\\慣性項\end{array}
\begin{array}{r}+~~~\underline{C\ v}~~~\\減衰項\end{array}
\begin{array}{r}+~~~\underline{K\ x}~~~\\剛性項\end{array}
\begin{array}{r}=~~~\underline{f}~~\\外力\end{array}
}$$
\begin{array}{l}
x&:& 変位(m) \\
v&:& 速度(m/s) \\
a&:& 加速度(m/s^2)\\
f&:& 外力(N)
\end{array}
\begin{array}{l}
M&:& 質量(kg) \\
C&:& 粘性減衰係数(Ns/m) \\
K&:& ばね定数(N/m)
\end{array}
Modelicaスクリプトの構成
上記の1次元のマスバネダンパーをModelicaで記述すると下記のようになります.
上側(緑色)が変数およびパラメータを宣言するエリア,下側(水色)が方程式を記述するエリアです.順に設定していきます.
スクリプト作成の手順
最初にModelicaスクリプトを書くための環境を用意するために,OepnModelicaを起動します.
モデルの新規作成
OpenModelicaを起動したら「メニュー>ファイル>Modelicaクラスの新規作成」もしくは をクリックしてモデルの新規作成をします.ここではモデル名をMyScriptとしました.左上のテキストビューのタブ を押してスクリプトを表示させると,下記のような「model (モデル名)」文と「end(モデル名);」文で囲まれたテキストが作成されますので,この2行の間にスクリプトを書いていきます.
model MyScript end MyScript;
ちなみに実際はannotationという文が挿入されていると思いますが気にしないでOKです(消してもいいです).単なる注記(アノテーション)ですので.
モデル新規作成の操作の詳細はこちらから
変数の定義
最初に変数を定義します.
\begin{array}{l}
x&:& 変位(m) \\
v&:& 速度(m/s) \\
a&:& 加速度(m/s^2)
\end{array}
下記のように実数を示すRealの後に変数を記載して宣言します.変数は英数字(記号はアンダーバーのみ可)であればもちろん1文字じゃなくても構いませんが,先頭の文字はアルファベットである必要があります.大文字小文字は区別されます.
Real x; Real v; Real a;
今回は単位は指定しません.自分の頭の中で把握しておくことにします.(単位の設定の仕方は後日アップ予定)
パラメータの定義
次に機械特性を定義します.
\begin{array}{l}
m&:& 質量(kg) \\
c&:& 粘性減衰係数(Ns/m) \\
k&:& ばね定数(N/m)
\end{array}
機械特性は通常パラメータですので,下記のように実数を表すRealの前に,パラメータである事を示すparameterを付けて宣言します.パラメータはシミュレーション中に動的に変化させることは出来ませんが,実行する度に数値の変更が可能です.
parameter Real m=1.0; parameter Real c=1.0; parameter Real k=1.0;
今回はシンプルにしたかったので数値はすべて1.0を入れました.(単位設定の仕方は変数と同様に,後日アップします)
方程式の定義
次に支配方程式を定義します.今回はシンプルに外力なしの自由減衰振動として$f=0$ とすると,支配方程式は $mx+cv+kx=0$ になります.速度は変位を微分したもの $v=\frac{dx}{dt}$ なので,1階の時間微分を表す関数 der を用いてv=der(x)と書きます.加速度は速度の微分 $a=\frac{dv}{dt}$ なので,下記のように3式を定義します.
equation v=der(x); a=der(v); m*a + c*v + k*x = 0;
ところで加速度は変位の2階微分でもあるので$a=\frac{d^2x}{dt^2}$とも書けますが,Modelicaには2階微分を表す関数はありません.
通常,他のほとんどの数値(数式)計算ソフト(Matlab,Mathematica,Maxima等)のシンボリック計算では高階の微分演算子が用意されてるので,なぜModelicaにはないのか疑問に感じるかもしれませんが,これは工学問題の支配方程式は機械系,電気系,流体,伝熱に関わらず,すべて2階以下であること,また高次の微分方程式であっても現代制御等で扱う状態空間表現にすれば,すべて1階微分の状態量で表現できること,などを考えると物理現象を記述するには1階微分の関数があれば十分です.
またModelicaは,モデル内のすべての連立方程式をすべて展開した後にコンパイルするので,1階微分で構成されている方がソルバーも実装しやすいでしょう.
モデルのチェック
すべてを記述し終わったら,メニュー>シミュレーション>モデルチェックもしくはアイコンをクリックしてスクリプトのチェックをします.
チェックが終わると下記のようなダイアログボックスが表示されます.
1行目は文法のチェックです.ここでsuccessと表示されずエラーが出た場合は,単純な文法ミスをしています.
2行目が最も重要で,方程式と変数の個数が表示されています.連立方程式を解くには変数の数と方程式の数が一致している必要があります.今回は変数と方程式は共に3個なので「3 equations and 3 variables」と表示されている事を確認します.
3行目は連立方程式における自明の式の数です.自明の式とは,その式だけで自ら明らかな(わざわざ方程式を解くまでもない)式のことです.このモデルでは v=der(x)と a=der(v) の2式は自明なので「2 of these are trivial equation(s).」と表示されています.
実行
モデルのチェックが終わったら,メニュー>シミュレーション>シミュレーションのセットアップ,またはアイコンをクリックして,終了時刻を30秒に設定してOKボタンを押します.
スクリプトが展開されコンパイルが始まります.しばらく待つとプロットモードに切り替わります.
結果を確認
プロット表示
プロットモードになったら,右側に「変数ブラウザ」が表示されますので,x にチェックを入れて変位を表示させてみましょう.すると・・・
変位が0と表示されています.またメッセージブラウザに
Assuming fixed start value for the following 2 variables: x:VARIABLE() type: Real v:VARIABLE() type: Real
というワーニングが表示されています.
このワーニングは「下記の2変数(xとv)の初期値に(適当な)推定値を入れました」という意味です.初期値を指定しなかったために表示された警告文です.計算としてはエラーではないのですが,初期値にデフォルトのゼロが入力されたため結果もゼロになっています.
初期値を指定
初期値を指定するために下記の文を追加します.これは,ばねの端部を(指やマニピュレータなどで)初期変位 x=1引っ張って離す動作を表しています.速度は与えない(静かに離す)ので初速度はv=0 です.
initial equation x=1.0; v=0.0;
この時この支配方程式は2階の微分方程式なので,3個の変数のうち2つに初期値を指定します.もしそれ以上(例えば加速度にも)与えると冗長になります.数値が適切であれば問題なく実行終了しますが,矛盾した値を入れるとエラーになるので注意してください.
もしくは下記のように,変数宣言と同時に属性として( )の中で与えても構いません.このようにModelicaの初期条件の与え方は二通りあります.
Real x(start=1.0); Real v(start=0.0); Real a;
ちなみに実務モデルでは,初期値は本当に必要な時以外は,ユーザー側では指定せずソルバーに推定させる設定にしておく事も多いです.特に複雑なアセンブリモデルの場合などは数値計算上の蓄積誤差等がエラーの原因になりやすいからです.
プロット表示(再)
もう一度実行すると下記のように表示されました.自由減衰振動の結果が表示されています.
ちなみに上記の波形は減衰が強いので,もしこの波形があまり「減衰振動」という気がしない・・・という方は粘性減衰係数をc=0.1程に下げてみてください.こちらの方が(いわゆる)自由減衰振動らしく感じられるでしょうか.
全スクリプト
最後に,今回のモデルの全スクリプトを下記に記します.
model MyScript Real x; Real v; Real a; parameter Real m = 1.0; parameter Real c = 1.0; parameter Real k = 1.0; equation v = der(x); a = der(v); m * a + c * v + k * x = 0; initial equation x=1.0; v=0.0; end MyScript;
今回はModelica(非因果的モデル)のエッセンスを紹介するために,できるだけシンプルなスクリプトにしました.
またModelicaはequation(方程式)エリアの記述において,連立方程式の順番や左辺と右辺を入れ替えても,実行結果に影響ありません(これを非因果的といいます).ただしパラメータを含む変数宣言のエリアは,方程式エリアであるequationより前におく必要があります.実際にいろいろ試してみると良いでしょう.