日本ソフトウェア科学会第 31 回大会 (2014 年度) 講演論文集

数式処理に基づくハイブリッドシステムシミュレータ

Hyrose の大規模モデルシミュレーションに向けた拡張
河野 文彦

小林 輝哉

松本 翔太

上田 和紀

ハイブリッドシステムとは時刻の経過に伴って状態や方程式系が離散変化したり連続変化したりするシステムである. HydLa は制約階層に基づいたハイブリッドシステム記述言語であり, Hyrose は数式処理に基づく HydLa の記
号実行シミュレータである. HydLa や Hyrose では多数のオブジェクトが相互作用するモデルを適切な効率で行うことが課題であった. そこで本研究では大規模なモデルを適切な効率で扱うために 2 つの拡張を施した. 1 つ目の拡張として, 大規模モデルの簡潔な記述を可能とするために集合記法の導入を行った. 2 つ目の拡張として, 制約と変数間の依存関係を利用して問題を複数の問題に分割して解き, 変数の連続性を利用して判定処理を省略することで大規模モデルに対するシミュレーション能力の向上を図った. 本論文では 2 つの拡張およびその評価について論じる.

であり, Hyrose では実行時間の増大が問題となって

1 はじめにいた.

ハイブリッドシステム [7] とは時刻の経過に伴って本研究では記述量の増大に対し, 集合記法の導入に状態や方程式系が離散変化したり連続変化したりす

よって対応した（4 節）. 実行時間の増大に対してはるシステムである.

次の 3 つの手法を用いることで対応した.

HydLa [11] は制約階層 [1] に基づいたハイブリッド

1. 制約と変数の関係グラフを構築し, 独立に計算システム記述言語であり, 制約間に優先度をつけるこ

可能な制約集合を特定し, 問題を独立な部分問題とでモデルの通常時の振る舞いと例外的な振る舞いとに分割する（5 節）.

を簡潔に記述することを可能としている. Hyrose [9]2. 上で求めた関係グラフを利用しつつ, 変数値のは数式処理に基づく HydLa の処理系であり, HydLa

連続性にも着目することで再計算が必要な制約で記述されたハイブリッドシステムの高信頼シミュを効率的に特定する（6 節）.

レーションを行う. 多数の自動車を対象とした交通

3. 関係グラフを利用して矛盾の原因となった制約モデルなど, 多くのオブジェクトが相互作用するモデ集合を特定し, その情報を利用し効率的に極大無ルの挙動は自明ではなく, そうしたモデルに対しての矛盾集合を求める（7 節）.

高信頼シミュレーションは有意義といえる. しかし

本研究では Hyrose で多数のオブジェクトが相互作

多くのオブジェクトが相互作用するモデルを HydLa,

用するモデルを適切な効率で扱うことを目的としてHyrose で扱う場合, HydLa では記述量の増大が問題いるため, オブジェクトの数の増大による実行時間への影響が確認しやすい 1 次元ビリヤード (図 1) のモ

Extensions for Large-Scale Simulation of Hybrid Systems Symbolic Simulator Hyrose
Fumihiko Kono, Teruya Kobayashi, Shota Matsumoto, 早稲田大学大学院基幹理工学研究科, Graduate School of Fundamental Science and Engineering,
Waseda University.
Kazunori Ueda, 早稲田大学理工学術院, Faculty of Science and Engineering, Waseda University.

デルを取り上げ, オブジェクトの増加に対応する手法の説明および評価を行う.

x1

図1

x2



xn

1 次元上での n 個の質点のビリヤード

2 ハ イ ブ リッド シ ス テ ム モ デ リ ン グ 言 語
HydLa
HydLa は制約階層に基づいたハイブリッドシステ
ム記述言語である. HydLa プログラム中の変数はす
べて時刻に関する関数変数となっており, HydLa プ
ログラムの意味は HydLa プログラムに記述された仕様を満たすような変数値の軌道である. HydLa では
制約を宣言し, 宣言された制約の優先度を宣言することによってモデルを記述する. HydLa の制約は比較演算子 ( “=”, “!=”, “<”, “<=”, “>”, “>=” ) を 1 つだ
け含む原子制約と論理記号, 時相演算子 “[]” からな

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

INIT(b,b0,vb0) <=> b=b0 /\ b’=vb0.
COL(b1,b2) <=> [](b1-=b2- =>b1’=b2’- /\ b2’=b1’-).
CONST(b) <=> [](b’’=0).
WALL(b) <=> [](b-=-1 \/ b-=10 => b’=-b’-).
INIT(x1,0,1),
INIT(x2,4,0),
INIT(x3,6,0),
INIT(x4,8,0),(CONST(x1), CONST(x2)) << COL(x1,x2),(CONST(x1), CONST(x3)) << COL(x1,x3),(CONST(x1), CONST(x4)) << COL(x1,x4),(CONST(x2), CONST(x3)) << COL(x2,x3),(CONST(x2), CONST(x4)) << COL(x2,x4),(CONST(x3), CONST(x4)) << COL(x3,x4),
CONST(x1) << WALL(x1),
CONST(x2) << WALL(x2),
CONST(x3) << WALL(x3),
CONST(x4) << WALL(x4).
図2

質点が 4 つの 1 次元ビリヤードの HydLa プログ
ラムる時相論理式によって記述される.
ド条件とする条件付き制約を表す. また変数の直後に

2. 1 HydLa プログラムの具体例“-” が付いたものはその変数の左極限値を表す. COL

本節では HydLa プログラムの具体例を示す. 図 1は 2 つの質点の位置が一致したときに質点同士が衝に示すモデルを HydLa で記述することを考える.

突し, 速度が入れ替わることを表している.

図 1 は n 個の質量の等しい質点が 1 次元空間上に4 行目では質点が等速運動をすることを表す制約並んだもので, 初期状態では質点 1 (x1) のみが正のCONST を定義し, 5 行目では質点と壁との衝突を表す

向きに動いており, 他の質点は静止している. また,

制約 WALL を定義している.

すべての質点の大きさは無視する. さらに, 質点同士

7 行目以降では定義した制約の宣言および制約の優の衝突および質点と壁との衝突の反発係数は 1 であ

先度の宣言をしている. HydLa では宣言した制約のるとする. 本モデルにおいて質点の数 n を 4 としたもことをモジュールと呼ぶ. 7 行目から 10 行目では各のを HydLa で記述した例が, 図 2 となる.

質点の初期値に関する制約を宣言している. 12 行目

1 行目は質点の初期値を表す制約 INIT の定義をから 17 行目では各質点の等速運動に関する制約およ行っている. 引数の b, b0, vb0 はそれぞれ質点の位び質点同士の衝突に関する制約の宣言および優先度

置を表す変数, 質点の初期位置, 質点の初速度を表しの宣言を行っている. この優先度の宣言により, 質点ている. 2 行目以降で使用している記号 “[]” は時相同士が衝突したとき, 衝突による速度変化が質点が等演算子 always であり, ∀t ≥0 を意味する. always が速運動をすることよりも優先され, 衝突した 2 質点のない制約は初期時刻のみで成立する制約である. 記号

速度が変化する. 19 行目から 22 行目では質点の等速“’” は時間微分を表すため INIT は質点の初期位置が運動に関する制約および壁と質点の衝突に関する制

b0 で初期速度が vb0 であることを表している.

約の宣言および優先度の宣言を行っている. この優先

2, 3 行目では 2 つの質点の衝突の制約 COL を定義

度の宣言では壁との衝突による質点の速度変化が質している. 引数の b1, b2 はそれぞれ異なる質点の位点が等速運動をすることよりも優先され, 壁と衝突し置を表す変数である. 記号 “=>” によって左辺をガーた質点の速度が変化する.

2. 2 解候補モジュール集合{1, 2, 3, 4}

宣言されたすべての優先度関係について式 (1) と{2, 3, 4}

式 (2) を満たすモジュール集合 MS のことを解候補{1, 3, 4}{1, 2, 4}{2, 3}{1, 4}{3}{2}{1, 2, 3}

モジュール集合と呼ぶ.
{3, 4}∀M1 ∀M2 (M1 <<M2 ∧ M1 ∈ MS
⇒ M2 ∈ MS)(1)∀R(¬∃M(R<<M) ⇒ R ∈ MS)(2){2, 4}{4}

HydLa では解候補モジュール集合のうち極大無矛盾{1, 2}{1}{}

集合が記述されたモデルの挙動を表す. 極大無矛盾集
合は時々刻々と変化する. 式 (2) は R よりも優先度が{1, 3}

図3

4 つの質点の 1 次元ビリヤードの解候補モジュール

高いモジュール M が存在しなければ R はすべての解集合の包含関係についてのハッセ図

候補モジュール集合の要素であることを意味し, このR のようなモジュールのことを required なモジュー
ルと呼ぶ.

• シミュレーションの記号実行を行うことができ,パラメータを持つシステムを扱うことが可能である

2. 3 解候補モジュール集合の具体例
図 2 のプログラムにおける全解候補モジュール集合
を集合の包含関係によるハッセ図で示すと図 3 とな

• パラメータの値によってシステムの振る舞いが定性的に変化する場合, 自動場合分けをしつつ網羅的に振る舞いの探索を行う

る. 図 3 では required なモジュールである INIT, COL,

Hyrose はシステムの離散変化時の振舞いを求める

WALL の記述を省略し, CONST(xn) を n と表記していポイントフェーズ (PP) と, 連続変化時の振舞いを求る. 以降, 要素がすべて数である集合はこの表記法にめるインターバルフェーズ (IP) とを繰り返し, 実行

よって表記されているものとする. 例えば {1, 2} は時に指定される終了条件（経過時刻やフェーズ数）が{INIT(x1,0,1), INIT(x2,2,1),満たされるまでシミュレーションを続ける. 各フェー

INIT(x3,4,1), INIT(x4,6,1),

ズにおいて, 各変数値が満たすべき制約の集合を無矛

COL(x1,x2), COL(x1,x3), COL(x1,x4),

盾性判定とガード条件判定により計算する. ガード条

COL(x2,x3), COL(x2,x4), COL(x3,x4),

件判定は条件つき制約に対するガード条件の成否判定

WALL(x1), WALL(x2), WALL(x3), WALL(x4),である. IP ではさらに次の PP の時刻の計算を行っ

CONST(x1), CONST(x2)}を表す.

ている．この計算では, 各ガード条件に対して成否が切り替わる時刻を計算し, その中で最小のものを次のPP の時刻として採用する. 最終的な Hyrose の出力

3 ハイブリッドシステムシミュレータ Hyrose

は, 各フェーズで求めた変数値の振舞いを時系列順にHyrose は HydLa の処理系であり, HydLa プログ

並べた軌道群となる.

ラムによって記述されたハイブリッドシステムの高信
頼シミュレーションを行う. Hyrose の特徴として, 以

4 集合記法の導入

下が挙げられる.

HydLa で多数のオブジェクトを含むモデルを記述• 数式処理を用いることで, シミュレーションの結果から誤差を排除している. 数式処理を行うために, バックエンドとして Mathematica [12] やする場合, すべてのオブジェクトに関する制約を列挙しなければならない.

1 次元ビリヤードの質点の数を n 個とした場合,

REDUCE [4] といった数式処理システムを利用

HydLa では質点同士の衝突の制約 COL の宣言およするび等速運動の制約 CONST との優先度の宣言は質点の全組み合わせ数である n C2 個の記述をしなければな
らず, 初期値の制約 INIT の宣言と壁との衝突の制約

WALL の宣言および CONST との優先度の宣言の数は n– 例えば図 4 の 1 行目の x1..x4 は {x1, x2,
x3, x4} を表す• Set1 or Set2
– 2 つの集合の和集合を表すとなる.
例として図 2 のプログラムでは WALL と CONST の宣言および優先度の宣言と INIT の宣言の数は 4 であ• Set1 and Set2
– 2 つの集合の積集合を表す

り, COL の宣言および COL と CONST の優先度の宣言

以上が集合の定義に用いる記法である. 次に HydLaの数が 4 C2 となっている.

プログラム内で集合を参照する際に用いる記法を以このように, オブジェクトの数が多い問題を HydLa
で記述しようとすると, 記述量が非常に多くなってし
まう. そこで HydLa に集合記法を取り入れることで記述量の増加に対応した.

下に示す.

• Name <=> Set– 制約の宣言および制約の優先度の宣言を要素とする Set を Name という名前で定義する

4. 1 導入した記法
多数のオブジェクトを含むモデルを簡潔に記述するために集合記法を導入する. 導入した記法において集合は数式の集合であるか, 宣言 (制約の宣言および制• Name := Set– 数式を要素とする Set を Name という名前
で定義する• Set[n]約の優先度の宣言の両方を指す) の集合である. 新た– 有限集合 Set の各要素に対して処理系が付に導入する記法では次のものを集合として取り扱う.

けた順番に基づき, その n 番目の要素を表す• {E1 , E2 , · · · , En }– E1 から En を要素とする集合 {E1 , E2 , · · · ,En } を表す– 各 Ei は数式で統一されているか, 宣言で統一されていなければならない• {E | SC1 , SC2 , · · · , SCn }– 条件式 SC1 から SCn をすべて満たすような E を要素とする集合を表す– 変数 Var が集合 Set の要素であることを意味する Var in Set と, 2 つの変数 Var1 と• |Set|– Set の要素数を表す• sum(Set)– Set の全要素の和を表す
さらに要素が制約の宣言および制約の優先度の宣言からなる集合を記述した場合, それは各要素の宣言であると解釈される. 例えば図 4 の 11 行目に書かれた集
合は {INIT(x[2],2*2-2,0), INIT(x[3],2*3-2,0),

INIT(x[4],2*4-2,0)} を意味しており, その集合の各要素を宣言したものであると解釈される.

Var2 が異なるものであることを表す Var1=!= Var2 の 2 つの式を条件式として使用で

4. 2 具体例

きる

導入した記法を用いた HydLa プログラムの例を図• s..e– 整数 s, e について {s, s + 1, · · · , e − 1, e}を表す– 例えば 1..3 は {1, 2, 3} を表す• Ns..Ne

4 に示す.
図 4 は図 2 と等価なプログラムである. 図 4 では

1 行目で各質点の位置を表す変数の集合を定義している. 3 行目から 8 行目は壁の位置を集合 x の要素数

|x| (この例では 4 となる) を用いて表していること以– アルファベット, 数字, “ ” からなる文字列

外は図 2 と同じである. 10 行目では集合 x の最初のN と整数 s, e を用いて {Ns, N (s + 1), · · · ,要素についての初期値の制約を宣言している. 11 行

N (e − 1), Ne} を表す

目は制約の宣言を要素とする集合の宣言を行っておれる. Modules(HP) は入力プログラム HP 中の全モ
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

x := x1..x4.
INIT(b,b0,vb0) <=> b=b0 /\ b’=vb0.
COL(b1,b2) <=> [](b1-=b2=> b1’=b2’- /\ b2’=b1’-).
CONST(b) <=> [](b’’=0).
WALL(b) <=> [](b-=-1 \/ b-=2*|x|+2=> b’=-b’-).
INIT(x[1],0,1),{INIT(x[i],2*i,0) | i in 2..|x| }.
{( CONST(x[i]), CONST(x[j]) )<< COL(x[i],x[j])| i in 1..|x|-1, j in i+1..|x| }.
{CONST(b) << WALL(b) | b in x }.
図4

導入した記法を用いた質点が 4 つの 1 次元ビリヤード HydLa プログラム

ジュールで, AtomicConstraints(m) はモジュール mに含まれるガード条件部を除いた原子制約の集合である. ガード条件はグラフの構築段階では考慮せず,

5. 2 節で述べる制約ノードの無効化において扱うものとする.

V は 変 数 ノ ー ド の 集 合 Variables(HP) で あ る.
Variables(HP) は, HydLa プログラム中に出現する全変数および時間微分の集合を表す.

Ep と E は互いに素なエッジの集合である．Ep は制約が変数の左極限値のみに言及していることを表すエッジ（以降 prev エッジと呼ぶ）の集合であり,Ep = {(c, v) | c ∈ C, v ∈ V, v.referredPrevBy(c)} と定義される. v.referredPrevBy(c) は変数 v の左極限値のみが制約 c 中に出現し, 変数値そのものには言及

り, 集合 x の最初の要素以外の初期値の宣言を意味ししないことを表す．

ている. 13 行目から 15 行目では質点同士の衝突の制

E は 制 約 か ら 変 数 値 そ の も の へ の 言 及 で あ り,

約と等速運動の制約の優先度の宣言を行い, 17 行目でE = {(c, v) | c ∈ C, v ∈ V, v.referredBy(c)} と定は壁と質点の衝突の制約と等速運動の制約の優先度

義される. v.referredBy(c) は c が変数 v の値そのもの宣言を行っている.

のに言及することを表す．なお, HydLa においては

図 4 のプログラムにおいて質点の数を増やすため

時間微分に関する制約は暗に微分される変数値の連には 1 行目で定義する集合の要素数を増やせばよい

続性も含意するため, 例えば制約 x0 = 0 は x と x0 のため, 集合の記法を取り入れた場合には質点の数が増両方に言及していることになる.

えてもプログラムの記述量は増大しない.

図 5 に, 1 次元ビリヤードプログラムにおいて質点を 3 つにした場合の関係グラフを示す. 図 5 のうち,

5 制約と変数の関係グラフ

長方形型のノードが制約ノードで, 楕円型のノードが本節では多数のオブジェクトの相互作用を含むモデ

変数ノードである. エッジについては点線が Ep の要

ルのシミュレーションを目的とした最適化手法につい素を表しており, 実線が E の要素を表す.

て述べる. このようなモデルは, 各オブジェクトの相互作用が発生する瞬間以外, 独立な部分問題に分割し

5. 2 関係グラフの利用てシミュレーションを進めることができる. この分割

前節で構築した関係グラフ上で非連結な成分は互いを行うため, 与えられた HydLa プログラムから制約に変数を共有しておらず, 矛盾し得ないため, シミュと変数の関係グラフを構築し, シミュレーションにおレーション中の無矛盾性判定やガード条件判定を独立いて利用する.

に行うことができる.
無矛盾性判定において結果が矛盾であると判定され

5. 1 関係グラフの構築た場合, 各制約ノードにはその制約が属するモジュー

本手法で構築するグラフは, 制約ノードと変数ノールの情報が含まれているため, 矛盾したモジュール集

ドからなる無向 2 部グラフ G = (C, V, Ep , E) である.

合を特定することができる. 特定したモジュール集合

C は制約ノードの集合で, C = {(c, m) | m ∈は, 7 節で述べる探索候補の導出に利用することがで

Modules(HP), c ∈ AtomicConstraints(m)} と定義さきる. また, ガード条件判定を行う際にも, そのガー

そのため, 本手法ではシミュレーション中の計算状態x1’’=0 (CONST(x1))

x1’’

時的に無効なものとして扱う. これにより非連結な成x1=0 (INIT(x1,0,1))x1’=-x1’- (WALL(x1))に応じて, グラフ中の制約ノードや prev エッジを一分を増やすことを可能とする.

x1

まず制約ノードを無効にできる場合とは, 以下の 3
条件のいずれかを満たす場合である.

x1’=1 (INIT(x1,0,1))x1’=x2’- (COL(x1,x2))

x1’

1. 制約ノードの属するモジュールが計算対象のモ
ジュール集合の中に含まれていない（ビリヤードの例における衝突時の CONST 制約

x1’=x3’- (COL(x1,x3))

がこれに対応する）

2. 制約のガード条件が成立していない

x2’=x1’- (COL(x1,x2))x2’’=0 (CONST(x2))
x2’=0 (INIT(x2,2*1,0))x2’’（ビリヤードの例における等速運動中の COL 制

x2’

3. 現在時刻がシステムの初期時刻でなく, かつ過

約や WALL 制約がこれに対応する）

去に有効になった always 制約でない

x2’=-x2’- (WALL(x2))

x2（ビリヤードの例における初期時刻以外の INIT
制約がこれに対応する）

x2’=x3’- (COL(x2,x3))

一方で prev エッジに関しては, PP 中は無効として良x2=2*1 (INIT(x2,2*1,0))

い. これは PP の左極限値が, 直前の IP の振る舞いから決定する定数値として扱えるためである．これに対x3’=x1’- (COL(x1,x3))

し IP 中は, 変数値と左極限値は必ず等しくなければ

x3’=x2’- (COL(x2,x3))x3’

x3’=0 (INIT(x3,2*2,0))

ならないため, prev エッジを有効なものとして扱う.

5. 3 関係グラフ利用の具体例

x3’=-x3’- (WALL(x3))

x3x3’’=0 (CONST(x3))

x3’’

本手法を用いて図 5 のグラフからノードやエッジを削減したグラフを, 図 6 と図 7 に示す. 図 6 は質点 1x3=2*2 (INIT(x3,2*2,0))と質点 2 が衝突する時刻の PP において, 全モジュー
ルを含むモジュール集合を計算対象とした場合の関係グラフで, 各質点ごとに独立な成分が得られてい図5

質点が 3 つの 1 次元ビリヤードの関係グラフ

る. この PP 中には, {COL(x1,x2), CONST(x1)} と{COL(x1,x2), CONST(x2)} がそれぞれ矛盾すること
ド条件中に含まれる変数と連結する制約ノードのみになるが, 本グラフを用いた無矛盾性判定を行うこと

を判定に用いればよいため, オブジェクトの増加によで容易にこれらの矛盾を検出することができる.

る計算時間の増加を抑えることができる.

図 7 は IP 中に有効な関係グラフで, 各質点ごとに

しかし HydLa プログラムから静的に関係グラフを

CONST 制約のみが関係する制約として得られている.

構築するだけでは, 図 5 のように全ノードが連結されこの関係グラフを用いると, x3 は x1 や x2 と無関係でた状態になることがほとんどである（本グラフ上で常あることがわかる. そのため, 例えば COL(x1,x2) のに非連結な成分とは, どの原子制約中でも互いに関わガード条件を判定する際に, CONST(x1) と CONST(x2)

りを持たない成分であり, 現実的なオブジェクトの相のみが必要であるという判断が可能になり, 扱う制約互作用をモデリングした場合にはほぼ発生しない）.

の数が削減される. 質点の数をどれだけ多くした場合

も各ガード条件判定の根拠となる制約は 2 つに減らx3x3’’=0 (CONST(x3))

すことができる.

x3’

6 計算結果の再利用

x3’’

本節では, Hyrose の実行アルゴリズム [9] のスケー
ラビリティ向上手法の一つとして実装した計算結果のx2’=x1’- (COL(x1,x2))x2’x2’’=0 (CONST(x2))

x2

再利用について述べる.
従来の Hyrose では, 制約のフェーズ間の差分と制約-変数間の依存関係を利用することで, 離散変化す

x2’’

る変数を特定し, 不要な再計算を排除していた [6].
離散変化する変数の特定において, 5 節で示した依

x1’=x2’- (COL(x1,x2))

x1’x1’’=0 (CONST(x1))

x1

存関係グラフを利用するよう実装を変更することで,離散変化する変数の特定に必要な制約集合の探索回数を削減した. また, ガード条件判定をさらに削減すx1’’

るために, ガード条件に含まれる変数の連続性に着目することによる判定の省略を新たに実装した. 以下,

図6

質点が 3 つの 1 次元ビリヤードの関係グラフ（質点

この 2 つについて説明する.

1 と質点 2 の衝突時）

x3x3’’=0 (CONST(x3))

x3’

6. 1 依存関係グラフを用いた再計算が必要な制約
およびガード条件の特定

5 節で述べた制約と変数の依存関係グラフを利用して離散変化した変数を求めるアルゴリズムを図 8 にx3’’示す. 図 8 内の記号 4 は集合の対象差を表す. アル
ゴリズムの基本的な工程は文献 [6] に示された Find-x2x2’’=0 (CONST(x2))

x2’ChangingVariables 関数と同様であるが, フェーズ間の差分に含まれる変数と関係した変数の探索に依存
関係グラフを用いる点が異なる. 以前の実装では現在

x2’’有効な制約の集合と直前のフェーズで有効な制約の集合とで共通する制約全てに対し, 離散変化した変数を

x1x1’’=0 (CONST(x1))

x1’
x1’’

持っているかを調べていたが, 依存関係グラフを用いて制約差分に含まれる変数と関係する部分グラフのみを探索するようにした.
例 え ば, n 質 点 の 1 次 元 ビ リ ヤ ー ド で , あ る 2
質点 (それらの位置を表す変数をそれぞれ xa, xb

図7

質点が 3 つの 1 次元ビリヤードの関係グラフ（連続変化中）と す る) の 衝 突 が 起 き た 場 合 , 以 前 の 実 装 で は 各

CONST(xk )(1 ≤ k ≤ n) に対して離散変化した変数を持っているか調べる必要があったが, 依存関係グラフ
を用いることで, 探索対象となる制約は CONST(xa),

CONST(xb) および COL(xa, xb) の 3 つに削減された.
また, 無矛盾性判定の対象とする離散変化した変Require: 直前のフェーズで有効な制約の集合 Sprev , 依存関係グラフ R
Ensure: 離散変化する変数の集合
1: S := R.GetConstraintStore()
2: CVtmp := CollectVariables(S 4 Sprev )3: CV := φ
4: for all v ∈ CVtmp do
5:
CV := CV ∪ R.GetRelatedVariables(v)6: end for7: return CV
図8

FindChangingVariables

このように変数の連続性を利用し, ガード条件判定
を省略する処理を実装した. 具体的には, 2 回目以降の IP において, 離散変化の原因ではない各ガード条
件に関して, 判定を行う前にそのガード条件に含まれる変数の時間微分値が直前の PP で定義されているか否かを調べる. 時間微分値が定義されてない変数がなければ, IP におけるガード条件の成否は直前の PPにおける成否と同一であると見なし, 判定を省略する.
再び, n 質点の 1 次元ビリヤードを考える. ある PP

数を含む制約の収集処理 (CollectChangingTell 関数

で 2 質点 (xa, xb) の衝突が起きた時, COL(xa, xb) の[6]) とガード条件判定の対象とする条件付き制約の収ガード条件 xa− = xb− が離散変化の原因となり, xa’

集処理 (CollectChangingAsk 関数 [6]) についても同と xb’ に関して新たに制約が追加される. 次の IP に様の改良を行った.

おいて, 従来の Hyrose では xa, xb を含む 2n + 1 個のガード条件を判定していたが, 変数の連続性を利用

6. 2 変数の連続性を利用したガード条件判定の省略
従来の実装では, 離散変化する変数と関連する全てした判定の省略を行うことで, 離散変化の原因でない
ガード条件の判定が省略され, 判定するガード条件は1 個になる.

のガード条件の判定を行っていたが, 離散変化する変数と関連するガード条件であっても, 変数の連続性を考慮することで判定を省略出来る場合がある. 具体的

7 解候補モジュール集合の動的な導出による探索候補の削減には, ある PP において離散変化の原因ではないガー

Hyrose は各時刻における極大無矛盾集合を計算す

ド条件に含まれる全ての変数が, その PP の時刻におることで HydLa プログラムのシミュレーションを行

いて連続である場合に, 次の IP における導出判定がう. 解候補モジュール集合の数が非常に多い場合, 愚省略出来る. その理由を説明するために, ガード条件

直な探索では非常に多くの時間がかかってしまう. そに含まれる変数 x1 , x2 , ..., xn に関する n 次元空間上

こで Hyrose の計算アルゴリズム上で計算が必要との, 時刻 t(t ≥ 0) における各変数の値を表わす軌跡

なった解候補モジュール集合を動的に生成することにg(t) を考える. ガード条件がある PP における離散変よって解候補モジュール集合の数の増大に対応する.

化の原因ではないとき, その PP の時刻を T とするHyrose の計算アルゴリズム上で計算が必要となっと, g(T) は, ガード条件が成立する部分空間と成立した解候補モジュール集合の集合のことを探索候補と呼ない部分空間との境界上には存在しない. また, ガー

ぶ. 本節では計算が必要となった解候補モジュール集

ド条件に含まれる変数値が T において連続であれば,合を導出し, 探索候補を更新しつつ極大無矛盾集合を

g(t) も T において連続である. よって, 次の IP の初求める手法について説明する.

期時刻 T + 0 における軌跡上の点 g(T + 0) は, PP の極大無矛盾集合の算出についての関連研究として

時刻 T における軌跡上の点 g(T) と同一の部分空間

文献 [8] がある. 文献 [8] では集合の要素となる制約にに存在する. このことは, PP と IP の初期時刻の間で優先度を設けていない点が本研究と異なる.

ガード条件の成否が切り替わらないことを意味しており, 同一の IP 内でガード条件の成否が切り替わる

7. 1 解候補モジュール集合の数ことはないので, PP におけるガード条件判定結果をモデルのオブジェクトの数が増えた場合, 解候補モ再利用することで, 次の IP におけるガード条件判定

ジュール集合の数が爆発的に増大する. 1 次元ビリを省略することが出来る.

ヤードでは質点の数 n に対し, 解候補モジュール集合

の数が 2n となる (図 3) . 例えば図 4 のプログラムでの全解候補モジュール集合は図 3 のようになっており, その総数は 24 個である.

HydLa プログラムではオブジェクトが 1 つのプロ
グラムにおいて解候補モジュール集合が s 個あった場
合, オブジェクトを n 個に増やすとオブジェクト間で相互作用がない場合のプログラムでは解候補モジュー
ル集合の数が sn 個となり, 相互作用がある場合にはそれ以上の数になる.

7. 2 解候補モジュール集合の動的な導出手法
図 9 は探索候補を更新するアルゴリズムである. 無
矛盾性判定の結果, 解候補モジュール集合が矛盾したときに呼び出される. 探索候補に含まれる解候補モ

Require: 今までの探索候補 MSS, 矛盾したモジュール集
合 IS, 極大無矛盾集合とわかっている集合 MCSS
Ensure: 新しい探索候補 N
1: N := ∅2: for all MS ∈ MSS do
3:
if MS ⊇ IS then4:
for all I ∈ IS do
5:
if ¬∃W ∈IS (W <<I ∧W 6≡I ) then6:
S := MS \ { M | M << I ∨ M ≡I }7:
if ¬∃ MCS ∈MCSS (S ⊂ MCS ) then8:
N := N ∪ {S}9:end if10:
end if
11:
end for
12:
else
13:
N := N ∪ {MS}14:
end if
15: end for
16: return N
図9UpdateMSS

ジュール集合は次の 3 つの性質を満たす.
ジュール集合の表記法で表したものである.

1. その時点で矛盾するかどうか不明である

ここでは質点 1 と質点 2 が衝突する瞬間の極大無2. 無矛盾とわかっている集合に包含されない
3. ハッセ図において, 矛盾した集合と直接繋がっ

矛盾集合を図 9 のアルゴリズムを使って求める.
まずは質点 1 と質点 2 が衝突した瞬間に極大無矛盾ている
ただし, 探索候補の初期値は, 全モジュールを要素と集合となる集合について確認する. 衝突の瞬間に質点する解候補モジュール集合を唯一の要素とする集合との衝突を表す制約 COL(x1,x2) のガード条件が成立すする.

るため, その後件の制約によって質点 1 と質点 2 の速

図 9 内の記号 ≡ は式 (3) で定義される.

M1 ≡ M2⇔

度が変化する. これは質点 1, 質点 2 が等速というこ(M1 <<M2 ∧ M2 <<M1 ) ∨とを表す制約である CONST(x1), CONST(x2) のそれ(3)

ぞれと矛盾する. 図 4 のプログラムでは (CONST(x1),

式 (3) は 2 つのモジュール M1 , M2 が同じモジュー

CONST(x2)) << COL(x1,x2) と優先度が宣言されていルであるか, 優先度が循環していることを表している.

るため, 質点 1 と質点 2 が衝突する瞬間では解候補モ

M1 = M2

優先度が循環していた場合, 式 (4) が成立し, M1 とジュール集合の中で, 無矛盾である集合 {3, 4}, {4},M2 がある解候補モジュール集合の要素であるか否か{3}, {} の中で包含関係において極大である {3, 4}は常に一致する.

が極大無矛盾集合となる. 図 3 において {3, 4} を強M1 ∈ MS ⇔ M2 ∈ MS

(4)

調して表示し, 矛盾する集合および矛盾する集合に繋

required なモジュールは全解候補モジュール集合のがっているエッジを薄く表示したものが図 10 となる.

要素であることがわかっているため, 探索候補の更新

次に図 9 を用いて極大無矛盾集合を算出する動作においては考慮する必要がない. 従ってこのアルゴリを示す. Hyrose は解候補モジュール集合のうち要素

ズムには required なモジュールを取り除いた集合を数が多い集合から無矛盾性を判定する. 従って最初に与えるものとする.

探索候補とする解候補モジュール集合はすべてのモジュールを要素とする集合である (図 11, Step 1).

7. 3 導出アルゴリズムの動作例

Hyrose は 5. 2 節で述べた関係グラフを利用し, 解候

図 4 のプログラムを用いて図 9 の動作を示す. 以

補モジュール集合を独立な複数の集合に分割して無矛盾

下, 要素が数のみで表される集合は図 3 での解候補モ性判定を行う. {1, 2, 3, 4} について無矛盾性判定を行












図 11



質点 1 と質点 2 が衝突したときに生成される解候補モジュール集合のハッセ図















成した集合が追加される. 以上で処理が終了して Nが返され, 探索候補が N に更新される. ここまでで生














































成された解候補モジュール集合のハッセ図は図 11 の

Step 2 となる.












































次に更新された探索候補 N ともう 1 つの矛盾集合{CONST(x2)} を用いて再度 UpdateMSS を呼び出す.
















同様の処理が行われて {3, 4} が返され, これが新たな探索候補となる. このときの生成された全解候補モ


ジュール集合のハッセ図は図 11 の Step 3 となる.
図 10

質点 1(x1) と質点 2(x2) が衝突したときの 1 次元
ビリヤードのハッセ図求められた矛盾集合すべてに関して UpdateMSS の呼出しが終了したため, 探索候補の中で要素数が最大である {3, 4} に対して無矛盾性判定を行う. この集

うと分割された集合のうち {CONST(x1), COL(x1,x2)}合は無矛盾であるため極大無矛盾集合となる. 図 11と {CONST(x2), COL(x1,x2)} がそれぞれ矛盾する集の Step 3 を見ると生成された解候補モジュール集合

合であることがわかり, UpdateMSS が呼ばれる. 引数の数が抑えられていることがわかる.

MSS に渡される集合は {1, 2, 3, 4} から required なモ
ジュールを除いた集合を要素とする集合 {{CONST(x1),

CONST(x2), CONST(x3), CONST(x4)}} であり, IS に

7. 4 1 次元ビリヤードにおける解候補モジュール
集合削減効果の考察{CONST(x1)} と {CONST(x2)} が, MCSS には空集合

本論文で用いた 1 次元ビリヤードのモデルでは初が渡される. 今回は矛盾する集合が 2 つあるため,期状態で 1 つの質点のみが動いており, すべての衝突

UpdateMSS を 2 回呼び出すことになる. まずは IS

に関する反発係数が 1 であり, すべての質点の質量がとして {CONST(x1)} を渡した場合を考える. 1 行目

等しいため, 同時に動く質点は 1 つのみである. さら

で N を初期化し, 2 行目で MSS 内の要素を取り出す.

に, 質点に大きさがなく, すべての質点の初期位置が

MSS の要素は {CONST(x1), CONST(x2), CONST(x3),

十分離れているため 3 つ以上の質点が同時に衝突しCONST(x4)} のみであるため, MS はこの集合になる.

たり, 2 箇所以上で同時に質点の衝突が起こったりす

MS が IS を包含するため 4 行目に移り IS の唯一のることはない. 従って, 1 次元ビリヤードモデルの状要素である CONST(x1) が I となる. W << I となる

態は次の 3 種類に分類できる.

ようなモジュール W が存在しないため 6 行目に移り,

1. 全質点が等速運動している

集合 {CONST(x2), CONST(x3), CONST(x4)} が生成さ

2. 1 つの質点が壁と衝突する

れる. MCSS が空集合であるため, 8 行目で N に生

3. 2 つの質点が衝突する





1. ではすべてのモジュール間に矛盾は生じない. 最初に無矛盾性判定を行う集合 (すべてのモジュールを

ジュール集合は生成されず, 1 つの解候補モジュール集合のみが生成される.

2. では質点 i が壁と衝突した場合, {WALL(xi ),

要素とする集合) が無矛盾であるため新しい解候補モ



CONST(xi )} のみが矛盾する集合となる. 2 つの質点が衝突したときと同様にして極大無矛盾集合 {1, 2,

· · · , i − 1, i + 1, · · · , n} が得られ, そのときに生成される解候補モジュール集合の数は 2 となる.

3. では上記の通り 3 つの解候補モジュール集合が




図 12

1 次元ビリヤードにおける第 3 フェーズ以降の各処生成される. 以上をまとめると各時刻に生成される解理の実行時間候補モジュール集合の最大数は 3 となり, この数は質点の数を増やしても変化することはない.

と考えられるが, 現状ではこれらの処理は各フェーズのボトルネックとなっていない.

8 評価実験

図 13 に質点の数 n の増加に対する Hyrose のシ

１次元ビリヤードに関する評価実験を行った. 本節

ミュレーション全体の実行時間の変化を示す. 全体ので扱う HydLa プログラムは図 4 をもとに, 質点の数

実行時間は O(n2 ) になっており, これはガード条件のを 10 から 100 まで 10 個刻みに変化させたものとし,

判定時間によるものだと考えられる. 前述の通り第 3

実行フェーズ数は 100 に設定した. 実験環境は表 1 にフェーズ以降の計算量については改善しており, 時間

示す.

がかかっているのは最初の IP におけるガード条件判

図 12 に, 最初の 2 フェーズを除いた, 第 3 フェーズ定である. 図 14 にガード条件判定の実行時間の全体

以降のガード条件判定および無矛盾性判定の実行時間に対する割合を示す. 質点の数が増えるごとに, 最初を示す. 無矛盾性判定もガード条件判定も, O(n) よの 2 フェーズの割合が増えていることがわかる. 6 節

り少ない計算量が実現できていることがわかる. これに述べたガード条件判定の省略は最初の IP を対象とらの判定処理においては, 5 節で述べた関係グラフのしておらず, 最初は n C2 個の COL 制約に関する判定

利用によって計算対象となる制約ストアのサイズを定を行う必要があるためである. この判定の計算量を減数オーダーとしている. 更に, 判定を行う回数についらすには, フェーズ間の差分情報の利用とは異なる方ては, 6 節で述べたフェーズ間差分の利用によって定針で取り組む必要があるが, それに関しては今後の課

数オーダーとしている. その上で完全な定数オーダー

題とする.

が実現できていないのは, 内部で用いているデータ構造の参照操作などが定数オーダーになっていないため

9 まとめ
本研究では, HydLa に集合記法を取り入れること

表1

実験環境

でオブジェクトの数の増加に対するプログラムの記述

OS

Ubuntu 14.04 LTS

量の増大を抑えた．また, 制約と変数の関係グラフを

CPU

Intel(R) Xeon(R) 2.50GHz,
Quad-core × 2 (8CPUs)

構築し, それを利用して問題を独立な小問題に分割し,Memory8GbyteBackendMathematica 9.0

不要な再計算および不要な解候補モジュール集合の削減を行った．1 次元ビリヤードの例題に関しては, 第3 フェーズ以降の実行時間を線形オーダー以内の計算



田研究室 HydLa 班の皆様に感謝する. 本研究の一部

は科学研究費補助金 (基盤研究 (B)26280024) の助成  



を受けて行った.








図 13

1 次元ビリヤードにおけるシミュレーション全体の実行時間





 


 





























図 14

1 次元ビリヤードにおけるガード条件判定時間の割合

量に抑えることができ, システムの初期状態の計算に関しても 2 乗オーダーにすることができた.
オブジェクトの数の増加による HydLa プログラムの記述量の増大および Hyrose の実行時間の増大は多数のオブジェクトが存在するモデルでは常に生じる問題であり, 本研究でその問題への対応を行ったことで,今後大規模ハイブリッドシステムのシミュレーションを行う上での足がかりを作ることができた.
今後の課題としては, 本手法を用いても効率化することができない最初の 2 フェーズの処理の効率化や,大規模モデルをシミュレーションする上での数式の複雑化への対応が挙げられる. 特に数式の複雑化は多くのモデルにおいて生じる問題であるため, 多様な例題を扱うためには数式の複雑さへの対応を行う必要がある.
本研究を行うにあたり, 貴重な意見をいただいた上参 考 文 献[ 1 ] Borning, A., Freeman-Benson, B., Wilson, M. :Constraint Hierarchies, Lisp and Symbolic Computation, Vol.5, 1992, pp.221–268.
[ 2 ] Carloni, L., Passerone, R., Pinto, A. andSangiovanni-Vincentelli, A. L. : Languages andTools for Hybrid Systems Design, Foundations andTrends in Design Automation, Vol.1, No.1, 2006,
pp. 1–204.
[ 3 ] Gupta, V., Jagadeesan, R., Saraswat, V. and Bobrow, D.: Programming in Hybrid Constraint Languages, in Hybrid Systems II, LNCS 999, Springer,
1995, pp. 226–251.
[ 4 ] Hearn, A. C., REDUCE,
http://reduce-algebra.com/.
[ 5 ] Henzinger, T. : The Theory of Hybrid Automata, in Proc. LICS 96, IEEE Computer Society
Press, 1996, pp. 278–292.
[ 6 ] 小林輝哉, 河野文彦, 松本翔太, 上田和紀： フェーズ
間の制約差分情報および制約-変数間の依存関係を用いた HydLa 処理系の最適化, 人工知能学会第 28 回全国大会, 4C1-2, 2014.
[ 7 ] Lunze, J. : Handbook of Hybrid Systems Control: Theory, Tools, Applications, Cambridge University Press, 2009.
[ 8 ] Malouf, R. : Maximal Consistent Subsets, Computational Linguistics, Vol.33, No.2, 2007, pp. 153–
160.
[ 9 ] 松本翔太, 上田和紀： ハイブリッド制約言語 HydLaの記号実行シミュレータ Hyrose, コンピュータソフト
ウェア, Vol.30, No.4, 2013, pp.18–35.
[10] 上田和紀, 細部博史, 石井大輔 : ハイブリッド制約
言語 HydLa の宣言的意味論, コンピュータソフトウェ
ア, Vol.28, No.1, 2011, pp. 306–311.
[11] 上田和紀, 石井大輔, 細部博史: 制約概念に基づくハ
イブリッドシステムモデリング言語 HydLa, SSV2008(第 5 回システム検証の科学技術シンポジウム), 2008.
[12] Wolfram Research, Inc., Mathematica,
http://www.wolfram.co.jp/products/
mathematica/index.html.