1Dモデルと単位系（Modelica）

前回の記事でマスバネダンパ構造を用いて，スクリプトの書き方を紹介しました．その際に説明をシンプルにするために単位をつけずに説明しましたが，Modelicaは物理現象を記述するために開発された言語ですから，Modelicaの中で扱う変数やパラメータには通常，単位を指定します．

Modelicaで単位を定義するには大きく分けて2種類の方法があります．「変数にそのまま属性として与える方法」と「Typeクラスをつかう方法」です．

変数宣言時に属性として与える
単位ミスによる重大事故
Typeクラスによる定義
1. 自分で定義する方法
2. SIunitsライブラリを使う方法
追記
1. import文でまとめる
2. 表示用の単位設定

変数宣言時に属性として与える

Modelicaの単位の記述法は非常に簡潔です．というのもModelicaには単位変換のアーキテクチャが標準で実装されているため，単位を指定するには変数宣言する際に，その属性として（　）内に「unit = “単位記号”」と単位をそのまま記述するだけでよいです．

変数

たとえば前回の記事の変数に単位を設定するためには，下記のように記述します．

\begin{array}{l}
x&:& 変位(m) \\
v&:& 速度(m/s) \\
a&:& 加速度(m/s^2)
\end{array}

Real x(unit = "m");
Real v(unit = "m/s");
Real a(unit = "m/s2");

パラメータ

パラメータも変数と全く同様に記述します．機械特性は下記のようになります

\begin{array}{l}
M&:& 質量(kg) \\
C&:& 粘性減衰係数(Ns/m) \\
K&:& ばね定数(N/m)
\end{array}

parameter Real m(unit="kg") = 1.0;
parameter Real c(unit="N.s/m") = 1.0;
parameter Real k(unit="N/m") = 1.0;

マスバネダンパーモデル

以上から，前回のスクリプトに単位を設定したものを下記に記します．

model MyUnit1

  Real x(unit="m");
  Real v(unit="m/s");
  Real a(unit="m/s2");

  parameter Real m(unit="kg") = 1.0;
  parameter Real c(unit="N.s/m") = 1.0;
  parameter Real k(unit="N/m") = 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 MyUnit1;

ほぼいつも通りに，そのまま書けばいいので簡単ですね．

プロット表示させてみると，今回は単位が表示されているのが分かります．

変数が単位表示されている

シンプルなモデルの場合はこれで十分です．しかしながら単位というのは，実機でのモデリングにおいてはもっと奥が深いものです．

単位ミスによる重大事故

多くの産業事故が単純なヒューマンエラーによって発生しています．

歴史上有名な単位計算ミスによる事故に，アメリカの火星探査機（Mars Climate Orbiter）の事故があります．1960年代頃（ケネディの時代）アメリカは多額の予算を投入して旧ソ連と熾烈な宇宙開発競争（Space Race）に明け暮れましたが，90年代になるともっと現実的であるべきとして，NASAでの宇宙開発もコストパフォーマンスや開発スピードを求められるようになってきました．（一般的な企業では当たり前ですが）

このような背景の元に起こったのが1999年の事故で，この事故は開発チーム間の単位系のやり取りのミスが原因と言われています．

Wikipedia（Mars Climate Orbiter）より

When NASA Lost a Spacecraft Due to a Metric Math Mistake | Simscale

失敗百選～単位系の取り間違いで火星探査機が行方不明～

「火星探査機事故の原因は測定単位系の混乱｜WIRED ニュース・レポート

この事故の損害額は約9400万ドルから1億2500万ドルとも言われており，このときNASAは大きく信頼を失いましたが，実は単位間違いというのは，シミュレーション・モデルの構築において，もっとも起きがちな深刻な凡ミスの一つです．

どんなに機械が高度に発達しても，その周辺にいるのは生身の人間ですから，ヒューマンエラーを一枚岩で防ぐことはできません．機械設計においてはこういったポカミスや危険源を現実的に回避するために，フェールセーフやフールプルーフなどのセーフティ対策を複層で行っています．

Modelicaにおいても，こういった単位ミスによるトラブルを防ぐための処方が備わっていますので，その使い方を説明していきます．

Typeクラスによる定義

自分で定義する方法

例えばモデルにダンパー機器を追加していくとします．ダンパーの粘性減衰係数をそれぞれ設定する事になると思いますが，ダンパー機器が増える度に毎回単位を記述するのでしょうか・・？

parameter Real c1(unit="N.s/m") = 0.1;
parameter Real c2(unit="N.s/m") = 0.2;
parameter Real c3(unit="N.s/m") = 0.3;
　　　　　・
　　　　　・
　　　　　・

簡単なモデルであればよいですが，複雑なモデルになると，どんどん変数やパラメータの個数は増えてくるので，それぞれに毎回記述していれば間違いの元です．よってModelicaでは変数やパラメータを指定する際には通常Typeクラスを使い，その中で単位定義を与えます．

※ クラスとはオブジェクト指向の概念の一つで，現代のプログラミングにおいては重要なものです．Modelicaはシンプルなオブジェクト構造を採用することで，システマティックなスクリプトを構成しやすくなっています．（オブジェクト指向についてはまた後ほど記事をアップしたいと思います）

Typeクラスを使うと，上記のダンパーの粘性減衰係数は下記のように定義することができます．

type Damper = Real(unit = "N.s/m");

parameter Damper c1 = 0.1;
parameter Damper c2 = 0.2;
parameter Damper c3 = 0.3;

ここでは粘性減衰係数のために「Damper」という名前のTypeクラスを作成しました．そうすると単位は一度だけこのTypeクラスの中で定義すればよく，複数の粘性減衰係数（ここではc1, c2, c3）を指定するには，Typeクラス（Damper）を参照して数値のみ指定すればよいです．

それではマスバネダンパーモデルに，Typeクラスを用いて変数定義した場合のスクリプトを下記に示します．

model MyUnit2 "Typeクラスを使った場合"

  type Position = Real(unit = "m");
  type Velocity = Real(unit = "m/s");
  type Acceleration = Real(unit = "m/s2");
  Position x;
  Velocity v;
  Acceleration a;

  type Mass = Real(unit = "kg");
  type DampingConst = Real(unit = "N.s/m");
  type SpringConst = Real(unit = "N/m");
  parameter Mass m = 1.0;
  parameter DampingConst c = 1.0;
  parameter SpringConst 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 MyUnit2;

ライブラリブラウザをみると，ツリー構造にTypeクラスが自動的に追加されています．

これでこのモデルは今後の拡張性が高くなりました．

しかし，ここでふと疑問が・・・単位って自分自身で定義せずとも，国際法で定められた規格があるのでは・・？

SIunitsライブラリを使う方法

そうです．Modelicaには国際規格の ISO/IEC 80000 に準拠したSI単位系のライブラリが存在しています．ライブラリブラウザの「Modelica＞SIunits」というライブラリを開いてください．

ここに約500ほどのSI単位がすでに定義されています．今回のマスバネダンパーモデルで使われている物理量のTypeクラスの記述を見ていきます．

変位

type Position = Length;
type Length = Real (final quantity="Length", final unit="m");

変位Position は長さのLengthと同じ次元なので，Lengthを参照するよう記述があります．

Type Lengthを見てみると，物理量（quantity）には長さ”Length”という属性が与えられ，単位（Unit）には “m” という属性が与えられています．また各quantityとunitの前のfinalという接頭辞は「変更を許可しない」という意味で，単位や物理量という根本的な属性がシミュレーション中に変更されないよう防いでいます．

速度

type Velocity = Real (final quantity="Velocity", final unit="m/s");

加速度

type Acceleration = Real (final quantity="Acceleration", final unit="m/s2");

質量

type Mass = Real (quantity="Mass", final unit="kg", min=0);

質量Mass には「min=0」とありますが，これは質量はマイナスの数値を取ることがないので最低値を0kgと宣言しているためです．もし間違ってマイナスの値を指定したりすると計算エラーになります．（粘性減衰係数やばね定数も，通常はマイナスの値をとる事はないのですが，特殊な材料開発の結果の可能性としてはありえます）

粘性減衰係数

type TranslationalDampingConstant=Real
       (final quantity="TranslationalDampingConstant", final unit="N.s/m");

ばね定数

type TranslationalSpringConstant=Real
       (final quantity="TranslationalSpringConstant", final unit="N/m");

これらのSIUnitsライブラリを使うと，スクリプトは下記のように記述できます．

model MyUnit3 "SIunitsライブラリを使用"

  Modelica.SIunits.Position x;
  Modelica.SIunits.Velocity v;
  Modelica.SIunits.Acceleration a;
 
  parameter Modelica.SIunits.Mass m = 1.0;
  parameter Modelica.SIunits.TranslationalDampingConstant c = 1.0;
  parameter Modelica.SIunits.TranslationalSpringConstant 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 MyUnit3;

単位は自分で定義せず，前定義された汎用的なSI単位系を使っているので，より安全度の高いスクリプトになっています．

追記

import文でまとめる

さらに通常Modelicaで変数を定義するときは，import文を使って簡略化する事が多いので，合わせて紹介しておきます．下記は Modelica.SIunitsライブラリをimport文でインポートしています．

model MyUnit4 "import文使用"

  import SI = Modelica.SIunits;

  SI.Position x;
  SI.Velocity v;
  SI.Acceleration a;

  parameter SI.Mass m = 1.0;
  parameter SI.TranslationalDampingConstant c = 1.0;
  parameter SI.TranslationalSpringConstant 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 MyUnit4;

Typeクラスを参照する際，毎回「Modelica.SIunits」と書いていましたが，「SI」の2文字で代用できるようになりました．

表示用の単位設定

ミスを防ぐためにも，計算実行時にはSI単位のMKS単位系（変位 m, 質量 kg, 時間 s）で統一すべきですが，通常慣例的に使っている単位がMKS単位系でないという事はよくあります．たとえば自動車の速度は通常 km/hですし，角度などは計算時はラジアン(rad)を使いますが，表示するときは一周360 °のディグリー(deg)であってほしいものです．

そういう際に最終的なプロット表示を，計算時とは別の単位で表示する方法を下記に示します．displayUnitという属性をつかうと，プロット時のデフォルト単位系を指定することができます．

下記は，速度vのプロットを cm/minで表示するように指示したスクリプトです．（Typeクラスの中で一括してdisplayUnit属性を定義することもできます）

model MyUnit5 "import文使用"

  import SI = Modelica.SIunits;

  SI.Position x;
  SI.Velocity v(displayUnit="cm/min");
  SI.Acceleration a;

  parameter SI.Mass m = 1.0;
  parameter SI.TranslationalDampingConstant c = 1.0;
  parameter SI.TranslationalSpringConstant 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 MyUnit5;

速度の表示が計算時と異なるcm/minで表示されました．変数ブラウザの単位をクリックすると「計算時の単位(m/s)」と「ディスプレイ用単位(cm/min)」の両方が選択できるようになっています．

表示用単位（cm/min）で表示されている

以上，シンプルな記述法から実務における使用法まで，単位の定義方法を説明しました．1DCAEの定義は「見通しのよい形式でシンプルに表現すること」です．要素開発時にざっと機能の確認をするとき，先行開発時にチームで作業するときなど，実機モデルのプロジェクトにおいては何を重要視するかで「見通しの良さ」は異なってきます．プロジェクト・リーダーは必要に応じた適切なモデルを構築していく必要があります．