システムとモデリング

SysML, Matlab/Simulink,他モデリング全般について。たまにプロジェクトマネジメント。

Modelicaプログラミング練習

こんばんは!Otepipiです! 今回はプログラミング言語Modelicaを使用してプログラミングの練習をしていきたいと思います。

Modelicaとは

Modelicaはオブジェクト指向プログラミング言語ですが、その中でも物理モデリングに特化しています。 言語仕様は下記にあります。

https://www.modelica.org/documents/ModelicaSpec34.pdf

大きな特長としては非因果モデルの記述が可能という点です。 非因果モデリングは、「出力」「入力」の区別をしない記述方法で、例えばJuliaやPython等で「3x+5=14+6α」といった式からxを求めたい場合

x = 3 + 2*α

とxについて式変形して求めることになりますが、Modelica言語ですと

3x + 5 = 14 + 6a

上記の式のままでxを求めることができます。

プログラミング練習

下記参考書から微分方程式の例題を説いてみます。

微分方程式で数学モデルを作ろう

微分方程式で数学モデルを作ろう

捕食ー被捕食者モデルを解いてみますが、これは以前に以下の記事でJulia言語で解いています。 otepipi.hatenablog.com

例題

$$捕食者-非捕食者の相互作用の単純なモデルは下記で与えられる。$$ $$\frac{dx}{dt}=x+y$$

$$\frac{dy}{dt}=-x+y$$

$$x………捕食者の個体数$$ $$y………非捕食者の個体数$$

$$初期値を下記のように仮定する。$$ $$x(0)=y(0)=1000$$

$$未来のxとyを決定せよ。$$

Modelicaでの実装

model TEST
  Real x (start = 1000);
  Real y (start = 1000);
 equation
 der (x) = x + y ;
 der(y) = -x + y ;
end TEST ;

さらにこのプログラムをコンパイルして実行するには下記の命令が必要です。

simulate(TEST, stopTime=1)

最後に結果をプロットします。

plot({x,y})

結果は下図のようになります。 f:id:Otepipi:20180221194825p:plain

参考:Juliaでの実装

参考に、Julia言語での実装を下に示します。

using ODE
using Plots
function f(t, v)
     (x, y) = v
     dx_dt=x + y
     dy_dt=-x + y
    
    [dx_dt ; dy_dt]
    end;

v_in = [1000.0; 1000.0];
time = 0.0:0.1:1;
x = map(v -> v[1], v);
y = map(v -> v[2], v);
plot(t,y,seriestype=:scatter,label="y")
plot!(t,x,seriestype=:scatter,label="x")

この時のプロットは下記のようになり、Modelicaの結果とよく一致しています。 f:id:Otepipi:20180221195501p:plain

今回は非因果モデリングの実例を示すことができなかったので 次回はModelicaを使った非因果モデリングについて解説したいと思います。

SysML モデリングソフトウェアを各種試したが定着しなかった話

こんにちは!Otepipiです。

前回の更新から2週間近く空いてしまいました。この間色々とSysMLのモデリングソフトウェアを試していたので纏めたいと思います。

前提条件

自分は今までこのブログに掲載していた各種SysMLダイアグラムをdraw.ioというウェブサービスを使って描画していました。

Flowchart Maker & Online Diagram Software

draw.ioはフリーで枚数制限無くダイアグラムが書けて、保存先も各種クラウドサービスやローカルを選ぶことができます。エクスポートできる拡張子も豊富にあるなど、フリーの描画ソフトでは群を抜いていると思います。 さらにSysMLの各種コンポーネントが予め用意されているのでダイアグラムの作成がとても捗ります。 なのでこのdraw.ioを特に不満なく使用していたのですが………

描画ソフトで作ったダイアグラムはモデルではない

 これまでdraw.ioでSysMLのブロック定義図やアクティビティ図を描いていましたが、これらはあくまで「図」で、「モデル」とは異なるものです。図を書くことは、真の意味のモデリングではありません。 自分がどれだけ理解しているかは怪しいですが、SysMLの言うところのモデルとは各種value、名前、関連などの要素で構成されたブロックまたはアクションそれ自体であり、動的なものです。  ブロック定義図やアクティビティ図はそれらモデル同士の関係を、ある観点から描画したものに過ぎず、絵である以上静的なものです。言うならばモデルの断面図と言った感じでしょうか?

そしてこのことからSysMLで言うところの「モデリング」作業には描画ソフトではなく専用のモデリングソフトウェアが必要になります。

SysML対応モデリングソフトウェア

各種SysML対応ソフトウェアをまとめていきます。

有料ソフトウェア

astah SysML

astah SysMLシステムモデリングツール | Astah

  • 日本企業によるソフトウェアで当然日本語対応
  • 1ラインセンス 12,000円/年
  • 評価版あり

Software ideas modeler

www.softwareideas.net

  • 日本語非対応
  • 1ライセンス 114ドル(永久ライセンス)
  • 評価板あり フリー版もあるにはあるのですが、これだとUMLしか使えません。 SysMLに対応するには有料のPREMIUMが必要になります。

Enterprise Architect

www.sparxsystems.jp

  • 日本語対応
  • 1ライセンス 79,400円(永久ライセンス)
  • 評価板あり

SysMLの使用にはシステムエンジニアリングverが必要で、それを購入する場合1ライセンス79,400円とかなり高額です。 個人向けでは無さそうですね。

無料ソフトウェア

Modelio

www.modelio.org

  • 日本語非対応
  • 無料

Papyrus

www.eclipse.org

この中から、特に無料ソフトウェアのModelioとPapyrusを導入してみました。

導入してみて

結論から言うと、ModelioもPapyrusも私のPC(win10 core i7 メモリ8GB)では動作がモッサリすぎてストレスが貯まる一方でした。 結局今でもdraw.ioのお世話になっています。私が作成しているようなモデルは複雑ではなく修正作業もほぼ無いので、「描画」で十分という面はあります。 モデリングソフトウェアは便利なのですが、動作が重いということが問題だと今回わかりました。 私のPCのスペックアップをすれば解決するかもしれません。

Julia言語で動的モデルのカオス挙動をシミュレーションする

こんばんは!Otepipiです。 今回も前回に引き続きJulia言語を使って数理モデリングとシミュレーションをしていきたいと思います。 例題は下記「数理モデリング入門」のものを使用しています。

数理モデリング入門 ―ファイブ・ステップ法― 原著第4版

数理モデリング入門 ―ファイブ・ステップ法― 原著第4版

シロナガスクジラの個体数xナガスクジラの個体数yの関係は以下で表すことができる。

$$ \frac{dx}{dt}=0.05x \left(1-\frac{x}{150,000}\right) -\alpha xy $$

$$ \frac{dy}{dt}=0.08y \left(1-\frac{y}{400,000}\right) -\alpha xy $$

初期値を下記のように仮定する。 $$x(0)=5,000 $$

$$y(0)=70,000$$

$$\alpha = 1.0×10^{-8}$$

未来のxyを決定せよ。

今回はPlotsパッケージを使用します。

using Plots

前回と同じ連立の常微分方程式ですが、今回はルンゲクッタ法ではなくオイラー法を使用します。下記は差分Δt=1年としたときのものです。 また、ステップ数はN=50回とします。

     const α = 1e-8
     const N = 50
     const Δt = 1
     x = zeros(N)
     y = zeros(N)
     t = zeros(N) 
     x[1] = 5000
     y[1] = 70000
     t[1] = 0
     
     for i = 1: N-1
        x[i+1]=x[i] + (0.05x[i]*(1-x[i]/150000)-α*x[i]*y[i])*Δt
        y[i+1]=y[i] + (0.08y[i]*(1-y[i]/400000)-α*x[i]*y[i])*Δt
        t[i+1] = i*Δt
       
                
            end


    
    

この結果をプロットします

plot(t,x,seriestype=:scatter,label="x", title = "Δt=1",xlabel="Year", ylabel="Population");
plot!(t,y,seriestype=:scatter,label="y");

f:id:Otepipi:20180206200400p:plain

シロナガスクジラの個体数xナガスクジラの個体数yともに解が現れています。このタイムスパンでは一定値への収束は確認できません。

次に時間間隔Δtを1年から10年に変更して再計算・結果をプロットしてみます。 f:id:Otepipi:20180206200425p:plain

解が一定に収束していくのを確認できます。 Δt=24ではどうでしょうか? f:id:Otepipi:20180206200437p:plain

ここでも解は収束するものの、序盤に振動が見られます。

では、Δt=27ではどうなるでしょうか? f:id:Otepipi:20180206200511p:plain

Δt=27になりますと、個体数は収束せず・発散もせず、周期2の振動を続けます。

Δt=32ではどうでしょう。 f:id:Otepipi:20180206200531p:plain

ここでも個体数は振動を続けますが、その周期は4に変わっています。

そしてΔt=37を見てみます。 f:id:Otepipi:20180206200558p:plain

ここではカオスが現れています。

Δtを37より大きくすると、解は無限遠に発散して計算できなくなります。

最後にこれらの画像をgifにしてみます。

まずΔt=1………37に対してオイラー法の計算を繰り返します。

     const α = 1e-8
     const N = 50
     const max = 37
    


     x = zeros(N , max)
     y = zeros(N , max)
     t = zeros(N , max)
     
for Δt = 1 : max
     x[1,Δt] = 5000
     y[1,Δt] = 70000
     t[1,Δt] = 0

     
     for i = 1: N-1
        x[i+1,Δt]=x[i,Δt] + (0.05x[i,Δt]*(1-x[i,Δt]/150000)-α*x[i,Δt]*y[i,Δt])*Δt
        y[i+1,Δt]=y[i,Δt] + (0.08y[i,Δt]*(1-y[i,Δt]/400000)-α*x[i,Δt]*y[i,Δt])*Δt
        t[i+1,Δt] = i*Δt
       
                
            end

end

    

次にΔt=1………37までのプロットをアニメーションにします。

anim = @animate for i = 1 : max

plot(t[:,i],x[:,i],seriestype=:scatter,label="x", title = string("Δt =", i),xlabel="Year", ylabel="Population")
plot!(t[:,i],y[:,i],seriestype=:scatter,label="y")

end
#gif(anim,"\anoi.gif",fps=30) エラーになってしまいました。

gifで書き出そうと思ったのですがgif()でエラーを吐いてしまいました。 今回は連番で書き出したpingをImageMagicで手動でGIF化することにしました。 完成したGIFがこちらになります。 f:id:Otepipi:20180206201232g:plain

解が収束した状態からカオスに至るまでが確認できました。

Julia言語を使って常微分方程式を解く(捕食者-非捕食者モデル)

こんにちは!Otepipiです。 今日はプログラミング言語Juliaを使用して微分方程式を解いてみます。 例題は下記「微分方程式で数学モデルを作ろう」から『捕食者ー非捕食者モデル』です。

微分方程式で数学モデルを作ろう

微分方程式で数学モデルを作ろう

例題

$$捕食者-非捕食者の相互作用の単純なモデルは下記で与えられる。$$ $$\frac{dx}{dt}=x+y$$

$$\frac{dy}{dt}=-x+y$$

$$x………捕食者の個体数$$ $$y………非捕食者の個体数$$

$$初期値を下記のように仮定する。$$ $$x(0)=y(0)=1000$$

$$未来のxとyを決定せよ。$$

Juliaでのプログラミング

まず、常微分方程式を解くためのODEパッケージとグラフ描画のパッケージを使用することを宣言します。

using ODE
using Plots

次に今回求める微分方程式を記述します。 コーティングの方法については下記サイトを参考にしています。

github.com

function f(t, v)
     (x, y) = v
     dx_dt=x + y
     dy_dt=-x + y
    
    [dx_dt ; dy_dt]
    end;

初期値とタイムスパンを指定します。

v_in = [1000.0; 1000.0];

time = 0.0:0.1:1;

常微分方程式を解きます。

t, v = ode45(f, v_in, time);

プロットの準備をします。

x = map(v -> v[1], v);
y = map(v -> v[2], v);

結果をプロットします。

plot(t,y,seriestype=:scatter,label="y")
plot!(t,x,seriestype=:scatter,label="x")

解析解との比較

結果が正しいか確認するため解析解と比較します。 この微分方程式系の解析解は以下になります。 $$ x = 1000e^t(\cos t + \sin t) $$ $$ y = 1000e^t(\cos t - \sin t) $$

今回のシミュレーション結果と解析解を比較します。

x_th=map(x -> 1000*exp(x)*(cos(x)+sin(x)), t);
y_th=map(x -> 1000*exp(x)*(cos(x)-sin(x)), t);
plot!(t,y_th,label="y_th")
plot!(t,x_th,label="x_th")

結果は良く一致していました。

カレーの作り方をSysMLでモデリングする(振る舞い図)

こんばんは!Otepipiです。本日は前回に続いてカレーの作り方をSysMLでモデリングしていきます。前回は構造図を作成しましたので今回は振る舞い図の作成となります。

 

前回の記事は↓

otepipi.hatenablog.com

カレーの作り方(再掲)

材料(4人分)

材料 分量
ハウス こくまろカレー<中辛>140g 1/2箱 70g
鶏肉(もも)   150g
玉ねぎ   中1個
じゃがいも   中1個
にんじん   中1/2本
サラダ油   大さじ1
  600ml

 

作り方

  1. 鶏肉、じゃがいも、玉ねぎ、にんじんは一口大に切る。
  2. 厚手の鍋にサラダ油を熱し、(1)の鶏肉、じゃがいも、玉ねぎ、にんじんをよく炒める。
  3. 水を加え、沸騰したらあくを取り、弱火~中火で材料が柔らかくなるまで約15分煮込む。
  4. いったん火を止め、ルウを割り入れて溶かし、再び弱火でとろみがつくまで約10分煮込む。

 

 アクティビティ図の作成

振る舞いを表現する図としてアクティビティ図を作成します。他SysMLの振る舞い図にはシーケンス図と状態機械図もありますが、食材の流れ(オブジェクトフロー)が重要な料理の工程ではアクティビティ図を使用するのが良いと思います。

f:id:Otepipi:20180201205724p:plain

ちょっとバランスの悪い図になってしまいましたが、上記カレーの作り方をアクティビティ図にすると上記のようになります。基本的にアクティビティ図は「左→右」もしくは「上→下」に流れるように作成します。またアクティビティ図の矢印は基本的に「物の流れを示す実線の矢印(オブジェクトフロー)」と「制御を示す破線の矢印(コントロールフロー)」の2種類がありますが、今回は食材という物の流れを示す工程なのでオブジェクトフローしか用いていません。また、図の枠線にかかっている四角(「鶏肉」「灰汁」など)は系外からの流入・流出を示します。

最後に前回掲載した構造図も示します。

 

f:id:Otepipi:20180128160906p:plain

料理では、構造図は「料理の材料」、振る舞い図は「調理方法」を示していると考えることができます。

このように普段何気なく使用している「レシピ」もモデリングすることができます。

 

カレーの作り方をSysMLでモデリングする(構造図)

こんにちは!Otepipiです。

モデリング能力を向上させるため、何か身近なものをモデリングしてきたいと考えています。今回はSysMLを使ってカレーの調理工程をモデリングしてみたいと思います。

レシピはハウス食品のWEBサイトにあったものを使用します。

housefoods.jp

 

材料(4人分)

材料 分量
ハウス こくまろカレー<中辛>140g 1/2箱 70g
鶏肉(もも)   150g
玉ねぎ   中1個
じゃがいも   中1個
にんじん   中1/2本
サラダ油   大さじ1
  600ml

 

作り方

  1. 鶏肉、じゃがいも、玉ねぎ、にんじんは一口大に切る。
  2. 厚手の鍋にサラダ油を熱し、(1)の鶏肉、じゃがいも、玉ねぎ、にんじんをよく炒める。
  3. 水を加え、沸騰したらあくを取り、弱火~中火で材料が柔らかくなるまで約15分煮込む。
  4. いったん火を止め、ルウを割り入れて溶かし、再び弱火でとろみがつくまで約10分煮込む。

 

 

構造図(BDD)

レシピの構成をブロック定義図で書き下したのが下記になります。

f:id:Otepipi:20180128160906p:plain

このブロックには食材しかありませんが、調理器具・食器・ガスや電気などのユーティリティーも本来であれば必要かもしれません(これはモデルの"粒度"の問題だと思いますが)。今回、各食材の量の指定が「~g」や「~個」、「大さじ~」など複数の表現に渡っていたのがモデリングの観点からすれば分かり辛い部分でした。この場合、例えば個数は『多重度』で指定することも可能ですが「~g」と言った表現は多重度で支持することが難しいと感じたため、多重度の使用はやめ、values で指定することにしています。

 

次回は振る舞いの表現になります。SysMLでは振る舞いの表現図としてアクティビティ図、シーケンス図、状態機械図の3つの図が用意されています。今回のような料理システムでは食材というオブジェクトがアクティビティに応じて変化していきますので、オブジェクトフローを扱うことができるアクティビティー図が適切だと思います。なので次回はアクティビティー図を用いた料理の振る舞いの様子をモデリングしていきたいと思います。

モデルをどのように記述するか?

こんにちは!Otepipiです!
前回は「モデルとは何か」について記事を書きましたので、今回は「モデルを記述する方法」について書きたいと思います。

 前回の記事は↓

otepipi.hatenablog.com

まずはモデルにはどのような種類があるのか見ていきたいと思います。
 

図形モデル

前回の記事でも登場しましたが↓は「宇宙」の図形モデルです。

 

f:id:Otepipi:20180114170254p:plain

このように実際の「宇宙」を抽象化して図形に落とし込むことができます。これが一般的な「モデル」のイメージでは無いでしょうか?最近は2次元の図形だけではなく3DCAD・3Dプリンタを使用して3次元図形でモデリングすることもできます。

 

 

数理モデル

↓は物質の拡散の様子を数式でモデリングした、有名なフィックの拡散法則です。


上記のように、「現象(フィックの拡散法則の場合は拡散現象)」のある一面を抽象化して数式で表現することができます。これは一般に日本語で「数理モデル」と呼ばれています。

 例えば小学校・中学校・高校のテストで出てきた所謂「数学(算数)の文章題」は、「日本語の文章を数式でモデリングする」作業を実質的にこなしています。数理モデルは我々にとっても馴染みのあるものです。

 

自然言語

自然言語とは日本語や英語など我々人類が使用している言語のことですが、自然言語によるモデルの記述とはつまり「文章」です。言葉による説明で事象をモデリングしています。先程数学の文章題の例を出しましたが

『ヒロシさんが家から6km離れたスーパーに時速3kmで~』のような文章題の「文章」自体が現象を自然言語で表現したモデルとなっています。

 

モデリング記述の問題点

通常、モデルはこれら3つ(図形・数式・自然言語)の複合となっていますが、これらのうち自然言語で記述されたモデルには大きな弱点があります。
・解釈の仕方がモデルを見る人によって異なる
ということです。文章でどのようにモデルを記述するか、については規定がなく各人によってバラバラですので、言葉が省略されていたり、違った意味で言葉が使われることが多くあります。

 

モデリング言語の登場

自然言語でのモデリングは解釈のバラつきが大きいため、バラつきの少ない言語でモデリングしたいというニーズがあり、そこで生まれたのがフローチャートUML、SysMLなどモデリング言語です。

 

今日はここまでに致します。