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

ハイブリッドシステム制約言語 HydLa の数式処理実行
系へのアフィン演算の導入
松本 翔太

上田 和紀

ハイブリッドシステムとは時間の経過に伴って状態が連続変化したり，状態や方程式系が離散変化したりする動的シ
ステムであり，物理学をはじめとした様々な分野への応用が可能である．HydLa はハイブリッドシステムの仕様を制約によって記述する宣言型言語であり，その処理系として Hyrose が公開されている．Hyrose は数式処理シミュ
レーションに記号定数を導入することで，パラメータを含むシステムの全振舞いを誤差なく求めることができるが，複雑な数式の出現による計算時間の爆発や，解析的に解くことのできない微分方程式を含むプログラムを扱えないことなど，数式処理に起因する問題を抱えている．精度を保証しつつこれらの問題を解決するために区間演算を用いる
ことが考えられるが，本研究では特に区間演算の拡張であるアフィン演算の手法を組み込むことを目標とする．ア
フィン演算は計算誤差の表現にダミー変数を使用することで，数値間の相関を考慮し，通常の区間演算で問題となる
区間幅の増加を抑制することを可能とする．ダミー変数は現在の Hyrose における記号定数として扱うことが可能であり，アフィン演算と Hyrose との親和性は高いと考えられる．本論文では，Hyrose にアフィン演算を組み込む際，新たに必要となる機能について論じるとともに，それらを用いたハイブリッドシステムのシミュレーションについて
論じる．

いる．

1 はじめにHyrose [9] †1 は HydLa プログラムのシミュレー

ハイブリッドシステム [8] とは時間の経過に伴ってタであり，ハイブリッドシステムの高信頼シミュレー

状態が連続変化したり，状態や方程式系自体が離散変ションを行うことで，その性質の理解や検証に役立て化したりするシステムを指す．ハイブリッドシステム

ることを目標としている．現実のシステムを適切にシのモデリングツールは物理学や制御工学，生命工学な

ミュレーションするためには，モデリングにともなうどさまざまな分野で重要な役割を果たすことができ

観測誤差や，システムの特性の環境や経年による変化

る．HydLa はハイブリッドシステムを制約によって等の不確定性を考慮する必要がある．従来の Hyrose

記述する宣言型言語であり，システムを表現する数式は数式処理に基づいた記号実行を行うことで，そのや論理式をそのままプログラムとして記述できるのような不確定性もパラメータとして扱うことを可能が特徴である．また，HydLa では制約間に優先順位にしていた．しかし記号実行を行う中で，パラメータ（階層）を設け，システムの通常時の動作と例外的なや無理数を含む数式はその長さの爆発的な増大や複動作の記述を簡潔に行うことを可能にしている．こ雑化を起こすことが少なくなく，Hyrose をより一般

れらの特徴により，HydLa ユーザが直感的かつ簡潔的なハイブリッドシステムのモデルに適用するためなモデリングを行うことが可能になると期待される．

には，これらの問題への対処が不可欠である．本研

HydLa の宣言的意味は文献 [14] において与えられて究では，誤差を厳密に扱う上でこれらの問題を解決するために区間演算 [12] およびその拡張であるアフィ

Shota Matsumoto, 早稲田大学大学院基幹理工学研究科
情報理工学専攻, Waseda University.
Kazunori Ueda, 早稲田大学理工学術院情報理工学科,
Waseda University.

ン演算 [7]（第 4. 2 節）の手法を Hyrose に組み込むこ†1 http://www.ueda.info.waseda.ac.jp/hydla/

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

2

14

INIT <=> 9 <= y /\ y <= 11 /\ y’ = 10.
FALL <=> [](y’’ = -10).
BOUNCE <=> [](y- = 15 => y’ = -4/5 * y’-).
INIT, (FALL << BOUNCE).

12

y

10

図1

天井に向かって投げ上げられる質点のモデリング例

86
4とを目標としている．アフィン演算は変数間の相関を2利用することで，計算結果の誤差の幅を少なく抑え

0
0.0

0.5

1.0

ることが可能な手法であり，アフィン演算で用いるダミー変数は現在の Hyrose においては記号定数としてそのまま扱うことができる．本稿ではアフィン演算を

1.5

2.0

2.5

3.0

t

図2

投げ上げられる質点の軌道

Hyrose に組み込む上で，必要もしくは実装すべき課

る舞いを表す BOUNCE（離散変化についての制約）を

題について考察を行い，更に Hyrose にアフィン演算

定義している．ここで，[]（時相演算子 always）はを組み込むための予備実験として，現在の Hyrose の

時刻 0 以降常に成り立つ制約を，=> は左辺を条件式

シミュレーションの途中で変数値をアフィン形式に近（ガード条件）とする条件つき制約を表す．最後に，4似した場合の影響について論じる．

行目で制約の半順序集合を構成している．
「,」は左右の制約の間に順序関係が存在しないことを，「<<」は

2 HydLa

右の制約が左の制約より優先されることを示す．こHydLa は制約を用いてハイブリッドシステムのモ

れらの記号でつながれた各要素が制約モジュールでデリングを行う宣言型言語である．システムが満たすある．この例では，FALL を BOUNCE より弱い優先度

べき個々の条件を，等式・不等式・微分方程式を原子で合成しているので，シミュレーション中は，FALL

論理式とする時相論理式で記述することでモデリン

より BOUNCE を優先して採用することになる．つまグを行う．HydLa プログラムは，このように時相論

り，FALL がこのシステムにおける通常時の動作であ理式で記述した制約（これを「制約モジュール」と呼り，BOUNCE が天井との衝突という例外的な動作であぶ）の集合である．各制約モジュールの間には優先順る，と解釈できる．

位を設けることができる．HydLa プログラム中に出こ の プ ロ グ ラ ム は INIT の 中 で ，y の 値 を現する変数は，時刻についての関数変数であり，それ9 <= y /\ y <= 11 と不等式制約で与えているの

らの変数の振る舞いのうち，プログラムの仕様を満た

で，y の初期値（以下 y0 とする）は 9 以上 11 以下

すものがプログラムの宣言的意味となる．

の任意の値が許され，解軌道は無限個考えられるこ

図 1 に，HydLa プログラムの例を示す．この例は，とになる．さらにこのモデルの場合，y0 が 10 以上か

高さ 15 の位置にある天井に向かって，その下から初否かによって BOUNCE 中の y- = 15 という条件がい

期速度 10 で投げ出される質点をモデリングしたもの

つか成立するか，もしくは永遠に成立しないかが変である．このプログラムに出現する変数は，質点の化する．また，y0 がちょうど 10 であるか否かによっ高さを表す y，およびその速度（一階微分），加速度て，離散変化時に FALL と BOUNCE とが矛盾するかど（二階微分）をあらわす y’，y’’ である．また，y-うかも変化する．このモデルの解軌道の一部をグラや y’- は y や y’ の左極限値を表す．1–3 行目は制約

フ化したものを，図 2 に示す（y0 を 9 から 11 まで，の定義であり，1 行目ではシステムの初期状態を表す

0.2 刻みで変化させた）．

INIT，2 行目では質点が重力によって下に引き寄せられることを表す FALL（連続変化についての制約），3行目では質点が天井に衝突する条件と衝突した際の振日本ソフトウェア科学会第 30 回大会 (2013 年度) 講演論文集

3

各計算の結果は，結果としてありうる実数の集合を包3 Hyrose

含する区間として定義され，例えば，加減算は以下のHyrose は HydLa の処理系であり，HydLa プログ

ようになる．

ラムによって記述されたハイブリッドシステムの高信

X + Y = [x + y, x + y]

頼シミュレーションを行う．現在の Hyrose の特徴と

X − Y = [x − y, x − y]して，以下が挙げられる．

区間演算は，計算機上での浮動小数点演算による計算• 一般的なシミュレータとは異なり，数値計算で誤差を考慮して，計算の結果を厳密に包含する用途はなく数式処理を用いることで，シミュレーショや，そもそも入力に誤差が含まれる場合の計算に用いンの結果から誤差を排除している．

られる．

• シミュレーションの記号実行を行うことができ，区間演算においては，計算によって得られた区間がパラメータを持つシステムを扱うことが可能で真の解を包含してはいるものの，その区間の幅が広がある．

りすぎることがある．こうした区間の拡大は，数学的• パラメータの値によってシステムの振る舞いがに同値であるはずの区間をそれぞれ独立なものとし定性的に変化する場合，自動場合分けに基づく探て扱うことで発生するものであり，区間の過剰な拡大

索を行う．

が結果の有意性を損なうことが少なくない．特にそのまた，Hyrose はバックエンドとして数式処理システ

ため，真の解を確実に包含しながら区間幅を小さく保ムである Mathematica [15] もしくは REDUCE [5] をつための計算手法の開発が求められており，次節で述用いるが，本稿では以降，Mathematica を利用する

べるアフィン演算はそうした手法の 1 つである．

場合に焦点を当てることとする．

Hyrose は離散変化時の振舞いを求めるポイント
フェーズと連続変化時の振舞いを求めるインターバ
ルフェーズ（以降それぞれ PP，IP と呼ぶ）を繰り返

4. 2 アフィン演算
アフィン演算は区間演算を拡張した演算であり，各変数の値をアフィン形式

し，実行時に指定される終了条件が満たされるまでシX = x0 + x1 ε1 + x2 ε2 + · · · + xn εnミュレーションを続ける．各フェーズで求めた変数値で表す．ここで εi はダミー変数と呼ばれ，−1 ≤ εi ≤ 1の振舞いをつなげた軌道群が，Hyrose の出力となる．

を満たす．各ダミー変数は独立な値をとるが，その影各 PP において値が 1 つに定まらない変数が出現し響の大きさは係数 xi が決定している．アフィン形式

た場合，Hyrose はその時刻での変数値に対応する記は以下のようにすることで通常の区間に変換するこ号定数を作成し，以降のフェーズではその記号定数とが可能である．

を参照して変数値の振舞いを求めていく．例として，[x0 − ∆, x0 + ∆], ∆ =

n
∑

|xi |

i=1

図 1 のプログラムに対しては，最初の PP で記号定例として，アフィン演算における加減算は以下のよう

数 y0 （ただし，9 ≤ y0 ≤ 11）を作成する．

に定義される．

4 区間演算とアフィン演算
本節では区間演算およびその拡張であるアフィン演
算について，概要とそれらの関係を説明する．

X ± Y = (x0 ± y0 ) + · · · + (xn ± yn )εnアフィン演算では計算対象となる各変数の値をア
フィン形式で保持することで，εi の各係数が計算過程で打ち消し合う可能性があるため，区間演算に比べて最終的な区間幅の増大を抑制することが可能と

4. 1 区間演算なる．一方でアフィン形式はダミー変数の線形結合に

区間演算は，実数もしくは以下のような実数の閉集なっているため，非線形演算を行う場合にはその線形

合 X （区間と呼ぶ）の間に定義された演算である．

近似を考える必要がある．その近似過程で生まれた誤X = [x, x] = {x ∈ R | x ≤ x ≤ x, x, x ∈ R}

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

4

差については新たなダミー変数を用意して対処する．

ことは可能である．近似を施すことで完全なアフィン演算を行う場合に比べて性能が劣る可能性はあるが，5 アフィン演算の導入に向けて近似を行う以前は変数間の相関を利用した計算が行

本節では Hyrose にアフィン演算を導入する際に，われており，また近似を行った後の計算も相関を利用

必要もしくは実装すべき事項についての考察を行う．

したものになる．そのため，常に相関を意識しない通常の区間演算と比較すれば，区間幅の増大を抑制でき

5. 1 アフィン形式の表現

ることが期待される．以下で，本稿の実験で用いた 2

まずアフィン演算を行うためには，各変数の値を第つの近似手法について説明する．

4. 2 節で述べたアフィン形式で表現できる必要がある．

5. 3. 1 0 次の近似

ここで，アフィン形式中のダミー変数 εi は，Hyrose

まず最も単純な近似の方法として，非線形項が出現において値の区間が [−1, 1] である記号定数（第 3 節）する数式を 1 つの区間にすることを考える．このたとして扱えばよいため，従来の Hyrose の実装を変更

めには，非線形項の中に出現するダミー変数を，通常することなく，アフィン形式を表現することが可能での区間演算における [−1, 1] の区間として考え，区間ある．

演算によって取りうる値の範囲を求めれば良い．この方法によって求められた区間は，ダミー変数が 1 つ

5. 2 非線形演算の近似のアフィン形式だと考えることができ，区間の中心をアフィン演算では変数値がアフィン形式を維持する

m，区間の幅を w としたとき，m + wεn+1 （n はこよう，非線形演算に対して線形近似を行って計算する

れまでに使用したダミー変数の個数）として扱うこと必要がある．一般の非線形単項演算や，一部の非線形になる．

二項演算に対しての最良の線形近似を行う方法 [11] は

5. 3. 2 1 次の近似

既に提案されており，それらの方法を Hyrose に組み

0 次の近似では元々含まれていたダミー変数を全て込むことができれば，Hyrose が内部的にアフィン演

削除してしまうため，その時点での変数間の相関が完算を行うことが可能になる．しかし，従来の Hyrose

全に失われてしまうことになる．これを避けるため，

はシミュレーションの中で，数式処理に基づいた非ダミー変数をできるだけ保ちつつ，非線形項を含む数

線形方程式の求解や，Quantiﬁer Elimination [1] を式を線形結合に近似することを考える．まず，元の数行っており，それらの中で行われている非線形演算を

式を f (ε1 , ε2 , · · · , εn ) とする．ここで f のマクロー

すべて線形近似された演算に置き換えるのは容易で

リン展開を考えると，ないと考えられる．一方で，一度でも非線形演算を許した場合，その結果には非線形な項が出現するため，
通常のアフィン演算の枠組みの中では扱いきれなく∂f (0, 0, · · · , 0)ε1
f (0, 0, · · · , 0) +∂ε1∂+··· +f (0, 0, · · · , 0)εn + R2∂εnのようになる（ただし，R2 はテイラーの定理におけなってしまう．そこで本稿では一時的に数式が非線形

る 2 次以降の剰余項とする）．R2 の大きさが無視できになることを許し，非線形な項を線形近似することでるほど小さいと仮定すると，上式から R2 を除いた式

この問題に対処することを試みる．その方法は次節でを f (ε1 , ε2 , · · · , εn ) の近似として扱うことができる．

説明する．

1 次の近似は 0 次の近似と比較して，ダミー変数を保持することで変数間の相関を有る程度残すことが5. 3 非線形項の近似できる反面，剰余項を無視する影響で近似結果が本来

非線形項が出現する数式をそのままアフィン演算のの結果を包含できていない可能性があるため，厳密に枠組みで扱うことはできない．しかし変数間の相関は誤差が扱えていないという問題点がある．この点にをある程度まで失うことを許せば，そうした数式をついては，剰余項 R2 を適切に評価し，その区間を包

通常の区間演算によって近似し，アフィン形式にする

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

5

含するためのダミー変数 εn+1 を使用することで解決できると考えられるが，それは今後の課題とする．

5. 4 アフィン演算を用いるべきかの判断基準
アフィン演算は通常の区間演算と比較して区間幅を抑えることができる手法ではあるが，発生する誤差は0 ではないので，数式処理で可能な範囲の計算は数式

INIT<=> x=3/4 /\ y=0 /\ vx=2 /\ vy=1/\ [](x’ = vx /\ y’ = vy).
RUN<=> [](vx’=0 /\ vy’=0).
BOUNCE <=> []((x-)^2+(y-)^2=1 =>vx = vx- - (x-*vx- + y-*vy-)*2*(x-) /\vy = vy- - (x-*vx- + y-*vy-)*2*(y-) ).
INIT, RUN<<BOUNCE.
図3

円の内側と衝突する質点の HydLa プログラム処理を用いた方が，最終的な計算結果の区間幅を抑えることができるはずである．このため，Hyrose におる場合を考える．数式処理シミュレーションによって求いてアフィン演算を用いるべきか否かの判断基準を適められた時刻 tg における x，y の値をそれぞれ xg ，yg

切に設定することが重要である．

とすると，x2g + yg2 = 1 は正確に成り立っている（x とそのような判断基準としては，数式処理を行う場

y はこの時刻でともに連続であるため，x− と x は同じ

合とアフィン演算を行う場合との計算時間の差や，ア意味を持つことに注意されたい）．しかし xg および yg

フィン演算を行った場合の結果の区間幅を反映していをアフィン形式で表すために区間値 Xa と区間値 Ya にるものが望ましいが，それらを事前に見積もるのは容近似したとすると，Xa2 + Ya2 ∋ x2g + yg2 = 1 は確実に易ではない．これまで Hyrose を用いて実験してきた

言えるが，同時に ∃xe ∈ Xa ∃ye ∈ Ya (x2e + ye2 > 1)

中でも数式がシミュレーション中にどのように複雑化や ∃xe ∈ Xa ∃ye ∈ Ya (x2e + ye2 < 1) も成り立ってし，またどういった場面で処理に時間がかかるのかをしまう可能性がある．このため，質点が円周を超えて見積もるのは，（問題の規模が大きい，といった単純な外に出ていく場合など，数式処理時点ではありえない

場合を除いて）困難であった．まず数式処理を試してと判断されていた状況に対しても，場合分けを行ってみて，あまりに時間がかかるようならアフィン演算にしまうことになる．このような事態を避けるためには切り替えて再実行する，という方法も考えられるが，ガード条件に直接関連している変数についてはアフィ

理想的には数式中の項の数や種類から，おおよその計

ン演算の対象としないか，もしくは近似を行った後で算時間を見積もることが可能となることが望まれる．

もそれ以前に成り立っていた制約の一部を，近似結果より優先されるべき制約として扱うことが必要とな

5. 5 数式処理によって得られた情報の区間演算における活用
数式処理によるシミュレーションとアフィン演算（および区間演算）によるシミュレーションとを切りる．前者の案によってアフィン演算を行えない変数が出現することは望ましくないため，後者の案を実装することが今後アフィン演算を組む込むに当たり特に重要な課題である．

替えて用いる場合，単にそれまで保持していた数式を区間に置き換えてシミュレーションを続行するだけで6 非線形項の近似を用いたシミュレーション

は，本来得られているはずの情報が失われてしまう可

本節では，従来の Hyrose のシミュレーションの中能性がある．
図 3 のプログラムを例とする．本プログラムは，質点が円の内側と衝突を繰り返すプログラムであで，第 5 節で述べた 2 つの近似手法を適用した場合の計算時間の変化や，解軌道の区間幅の変化について述べる．

り，質点と円との接触を判定するガード条件としてこれらの実験は，将来的に Hyrose に正式なアフィ(x-)^2 + (y-)^2 = 1 を含んでいる．ここで，Hy-ン演算を組み込むことに向けた予備実験として行ってrose の数式処理シミュレーションによって時刻 tg で前おり，数式処理シミュレーションの中で区間演算を用

述のガード条件が成り立つことが判明し，その時刻かいることで，結果の精度や計算時間に現れる影響を確らアフィン演算を用いたシミュレーションに切り替え

かめることを目的としている．実験環境は，OS: De-日本ソフトウェア科学会第 30 回大会 (2013 年度) 講演論文集

6

1フェーズあたりのフェーズあたりのフェーズあたりの時間
時間 ( μsec )<=> y = 10 + e1 /\ y’ = 0 + e2/\ e = 4/5 + e5/\ -10^(-5) < e1 < 10^(-5)/\ -10^(-5) < e2 < 10^(-5)/\ -10^(-5) < e3 < 10^(-5)/\ -10^(-5) < e4 < 10^(-5)/\ -10^(-5) < e5 < 10^(-5)/\ [](e’ = 0 /\ e1’ = 0 /\ e2’ = 0/\ e3’ = 0 /\ e4’ = 0 /\ e5’ = 0).
FALL<=> [](y’’ = -10 + e4).
BOUNCE <=> [](y- = 0 + e3 => y’ = -e * y’-).
INIT, FALL << BOUNCE.

1E+09

INIT

100000000

10000000
近似無し
0次近似

1000000

1次近似10000010000
PP IP PP IP PP IP PP IP PP IP PP IP PP IP PP IP PP IP PP IP
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

図5
図4

図 4 の実行時間比較

地面と衝突する質点のプログラム
0.00120.001

2.3GHz，メモリ: 16Gbyte，使用した Mathematica

0.0008
時刻
時刻のの区間幅

bian（etch），CPU: Quad- Core AMD Opteron(tm)のバージョン： 8.0 となっている．

近似無し
0.00060次近似
1次近似

0.0004

6. 1 地面と衝突する質点のプログラム
0.0002

1 つめの例題として，地面との衝突を繰り返す質点0のモデルに対し，質点の初期位置と初期速度，重力加速度，地面の高さ，衝突時の反発係数の 5 つの値にパ

PP3

PP5

図6

PP7

PP9

PP11 PP13 PP15 PP17 PP19 PP21

図 4 の実行結果の区間幅比較

ラメータを含めたプログラムを考える（それぞれ e1，e2，e3，e4，e5 とする）．本プログラムは，アフィン

験では PP11 以降のフェーズのシミュレーションは，演算におけるダミー変数を想定したパラメータを含時間がかかり過ぎるため打ち切っている．近似を行っ−5≤ εk ≤ 10

−5

としている

た場合については，0 次近似でも 1 次近似でも計算時（理由については，第 5. 3. 2 節で述べている）．各制

間に大差は無いが，近似無しの場合に比べて計算時間

約の意味は図 1 と同様である．

を抑え，シミュレーションを先に進めることができてんでおり，その変域を 10

本プログラムについて，近似無し，0 次近似あり，いることが分かる．

1 次近似ありのそれぞれの条件でシミュレーションを次に図 6 の区間幅の比較について考える．近似を行った．なお，近似を行った場合については最初の衝行った場合でも，そうでない場合と比べて区間幅の拡突に対応する PP3 において，y’ に関する近似を行っ大に大差は無いが，0 次近似を行った場合は 1 次近似ている．各フェーズの実行時間の比較を図 5 に，区間を行った場合に比べて若干区間幅が広くなっている．

幅の比較を図 6 に示す．図 6 では，各 PP の時刻と一方で 1 次近似を行った場合は，僅かだが近似無し

して得られた数式の意味する区間幅を比較している．

の場合よりも区間幅が狭くなっており，これは近似を

時刻の区間幅を比較するのではなく，質点の速度や質行ったために真の解を包含しきれなくなったことを点の位置について比較することも考えたが，そもそも意味している．この理由として，本稿で行っている 1の衝突時の時刻が異なってしまう状態での比較になっ次近似が 2 次以降の項の計算を打ち切っており，そてしまうため，本実験では時刻の区間幅を比較対象との誤差を正確に扱っていないという点が考えられる．

して採用した．

こうした誤差を扱うために，打ち切った項の誤差評価

図 5 を見ると，近似無しの場合は，近似を行った
場合と比べてフェーズが経過するごとに 1 フェーズ
あたりにかかる時間が指数的に増大している．本実を精確に行う必要があるということが分かった．

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

図7

4

2

0

y

INIT <=> x = 0 /\ y = 0 /\ x’ = 1 + e1/\ y’ = 2 /\ -10^(-12) <= e1 <= 10^(-12)/\ [](e1’ = 0).
CONST <=> [](x’’ = 0 /\ y’’ = 0).
BOUNCE(tan, d) <=>[](y -= tan * x- - d =>x’ = (-1 * tan^2/(1+tan^2) + 1/(1+tan^2))* x’- + (2 * tan/(1+tan^2)) * y’/\ y’ = (2 * tan/(1+tan^2)) * x’+ tan^2/(1+tan^2)+ -1/(1+tan^2)) * y’-).
INIT, CONST<< (BOUNCE(3^(1/2), 2), BOUNCE(1,-2),
BOUNCE(-3^(-1/2), 3), BOUNCE(-1, -4)).

7-2四角形の各辺と衝突する質点のプログラム

6. 2 四角形の各辺と衝突する質点

-4
-4

-2

0

2

4

x

2 つめの例題として，2 次元空間において質点が不等辺四角形の内部より等速度運動を開始し，四角

図8近似を行わない場合の図 7 の解軌道

形の各辺との衝突を繰り返すというモデルを扱う．

HydLa プログラムを図 7 に示す．質点の位置は水平4方向が変数 x，垂直方向が変数 y によって表されており，k が反発係数を表している．制約 INIT はシステ
ムの初期状態を，CONST は等速度運動を表しており，

2

BOUNCE は引数つき制約として定義されており，それる．本プログラムではパラメータ e1 を使用しており，y

ぞれの呼び出しが四角形の各辺との衝突を表してい010−12 ≤ e1 ≤ 10−12 としている．前節の例題よりパラメータの変域を小さく設定したのは，誤差の影響による定性的な軌道の変化を回避し，各手法の比較を定-2

量的に行いやすくするためである．
本例題を Hyrose でシミュレーションする際に近似を用いる場合，第 5. 5 節において述べたような理由か

-4
-4

-2

ら，質点の位置に関する近似を行った際に場合わけが発生してしまい，シミュレーションがかえって複雑化

0

2

4

x

図9

0 次近似を行った場合の図 7 の解軌道してしまう結果となった．そのため，本稿では第 5. 5
節の課題を解決することができたと仮定した場合のシ

実行時間の比較を図 11 に，区間幅の比較を図 12 にミュレーションへの影響を見積もるため，図 7 と同示す．本実験は Hyrose のシミュレーションと異なり，様の軌道を求めるためのプログラムを Mathematica

離散変化を行うフェーズと連続変化を行うフェーズとによって実装し，各近似の効果を計測した．本実験でを区別していないため，図 11 と図 12 の横軸は質点の近似は全衝突時刻において，全変数について行ってと辺との一回の衝突に対応していることに注意されいる．

たい．また，図 8 と図 10 ではパラメータの影響が微近似無し，0 次近似，1 次近似のそれぞれの場合の軌道を図 8，図 9，図 10 に示す．また，各手法での小であるため，軌道はほとんど一本の折れ線として描画されている．

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

1 ´ 105
5 ´ 104Byte

2

1 ´ 104

y

5000

1000

0

5

10

15

20

25

30

Phase

図 13

-2

図 7 の実行時の数式サイズの比較（赤：近似無し，青：0 次近似，緑：1 次近似）

-4
-4

0

-2

2

4

x

図 10

近似無しの場合は，10 回目の衝突までのシミュレー
ションができたが，それ以降のシミュレーションは実

1 次近似を行った場合の図 7 の解軌道

行時間が指数的に増大したため確認することができな
かった．このように実行時間が増大する理由として数式のサイズの増大が考えられる．本例題の実行時の数

100

Tim e H secL

式のサイズの変化を，図 13 に示す． 図 13 では，各フェーズにおける変数値（x, y, x′ , y ′ の 4 つ) を表す数

10

式のサイズを，Mathematica 組み込みの ByteCount1関数によって調べた．ByteCount 関数は，引数として与えた式を Mathematica 内部で記録するために必

0.1

要なバイト数を返す．図 13 を見ると実行時間と同様
5

10

15

20

25

30

Phase

図 11

に，近似を行った場合は数式のサイズがある範囲に収まっているのに対し，近似無しの場合は指数的に増大

図 7 の実行時間の比較（赤：近似無し，青：0 次近似，緑：1 次近似）

しているのが分かる．

0 次近似を行った場合は 35 回目の衝突までのシミュ
0.1

レーションが実行できたが，36 回目の衝突において区間幅の増大により場合分けが必要になった．実装しWidth of Interval

0.001

た Mathematica プログラムでは場合分けに対応して10-5

いないため，それ以降のシミュレーションは行っていない．34 回目や 35 回目の衝突では場合分けこそ発生

10-7

しないが，区間幅が拡大しすぎて有意な結果とは言い10-9難いため，本稿では 33 回目の衝突までの結果を掲載
5

10

15

20

25

30

Phase

図 12

図 7 の区間幅の比較（赤：近似無し，青：0 次近似，緑：1 次近似）

する．図 9 は各衝突時刻で質点の位置や速度について，下限から上限まで 0.05 刻みで動かしたものを描画している．質点が四角形の外に飛び出しているように見える部分があるが，これは近似の結果として質点の位置がずれたためであり，定性的には四角形の内部で跳ねる軌道として扱われている．0 次近似を行う

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

9

場合は近似を行うたびに指数的に区間幅が増大していってしまい，あまり長いシミュレーションは有意性

8 まとめと今後の課題を失ってしまうが，近似を行わなかった場合に比べて本稿では HydLa の記号実行シミュレータである

実行時間の増大を抑えることはできている．指数的にHyrose に対し，新たにアフィン演算を組み込むため区間幅が増大する理由としては，変数間の相関を考に必要もしくは実装すべき事項についての考察を行慮しない 0 次近似を繰り返し行っているためであり，い，組み込みに向けての一歩として非線形項の線形0 次近似を行う回数を減らすことで，ある程度は増大項への近似を提案した．更に，予備実験として現在のを緩和できると考えられる．

Hyrose のシミュレーションの中で線形項への近似を

1 次近似を行った場合には近似無しの場合や 0 次近行った場合の，実行時間や誤差の影響による区間幅の似の場合のような限界には達しない．8 回目の衝突ま変化を調べた．実験の結果として，シミュレーション

では 0 次近似の場合と比べて区間幅を抑制することの実行時間を抑制することはできたが，区間幅の増大ができており，かつ近似無しの場合よりも実行時間のを抑制しつつ真の解を完全に包含することはできて増大を抑制することができている．特に 8 回目の衝

おらず，そのためには 1 次近似において計算を打ち突においては，1 次近似を用いた場合は 0 次近似と切っている剰余項の評価を行うことが重要であること比較して区間幅を 1/1000 に抑え込むことができておが分かった．今後の課題としては，剰余項を精確に評り，変数の相関を利用することの効果が表れていると価し包含する方法や，ガード条件に出現する変数に言える．一方で 1 つめの例題と同じく，この問題にお対して近似を行っても，シミューレション上で無駄ないても 1 次近似を行った場合の軌道は，近似無しの場合分けを発生させることが無いような方法の提案，場合の軌道を完全には包含できておらず，その乖離は実装を行いたいと考えている．

特に 9 回目の衝突から顕著になっている．実行時間

謝辞 本研究を行うにあたり，貴重な意見を頂いた上の増大と区間幅の増大の抑制を両立するために，1 次田研究室 HydLa 班の皆様に感謝する．また，本研究

近似において計算を打ち切っている項の誤差を正しくの一部は, 科学研究費補助金 (基盤研究 (B) 23300011)

扱うことは重要だと考えられる．

の補助を得て行った.

7 関連研究本研究と同じく，区間演算を用いてハイブリッドシ
ステムのシミュレーションを行う研究として，Hybrid

cc [3] [4] の実装や Hyson [10] が挙げられる．Hybrid
cc では HydLa と同じく制約を用いた簡潔な記述がで
き，区間演算による精度保証も行われているが，離散変化を厳密に扱うことはできていない．Hyson は不確定性を持つ Simulink プログラムに対して区間演算を用いてシミュレーションを行うツールであるが，現時点では解軌道の精度保証は行われていない．これらの研究と比較すると，本研究は数式処理と区間演算を併用している点と，解軌道の厳密な精度保証を目指している点が異なっている．

10
参 考 文 献

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

stolﬁ.html#and-com-sto-94-aax[ 8 ] Lunze, J. : Handbook of Hybrid Systems Con[ 1 ] 穴井宏和，横山和弘：QE の計算アルゴリズムとそ
trol: Theory, Tools, Applications, Cambridge Uniの応用，東京大学出版会，2011.
versity Press, 2009.
[ 2 ] Carloni, L., Passerone, R., Pinto, A. and[ 9 ] 松本翔太，上田和紀：ハイブリッド制約言語 HydLa
Sangiovanni-Vincentelli, A. L. : Languages andの記号実行シミュレータ Hyrose，コンピュータソフト
Tools for Hybrid Systems Design, Foundations andウェア，採録決定済み
Trends in Design Automation, Vol. 1 No. 1, 2006,[10] Mimram, S., Bouissou, O. and Chapoutot, A. :pp. 1–204.
HySon: Set-based Simulation of Hybrid Systems, to[ 3 ] Carlson, B. and Gupta, V. : Hybrid cc with Inbe presented at IEEE International Symposium on
terval Constraints, in Proc. HSCC 98, LNCS 1386,
Rapid System Prototyping (RSP’2012), 2012.
Springer, 1998, pp. 80–94.
[11] 宮島信也，宮田孝富，柏木雅英：アフィン演算に[ 4 ] Gupta, V., Jagadeesan, R., Saraswat, V. and Boおける最良乗算について， 電子情報通信学会論文誌，
brow, D.: Programming in Hybrid Constraint LanVol.J86-A，No.2，2003, pp. 150–159，
guages, in Hybrid Systems II, LNCS 999, Springer,[12] Moore, R. E., Kearfott, R. B., Cloud. M. J.: In1995, pp. 226–251.
troduction to Interval Analysis, Society for Indus[ 5 ] Hearn, A. C., REDUCE,
trial and Applied Mathematics, 2009.
http://reduce-algebra.com/.
[13] Tabuada, P. : Veriﬁcation and Control of Hybrid[ 6 ] Henzinger, T. : The Theory of Hybrid AuSystems，Springer-Verlag, 2009.
tomata, in Proc. LICS 96, IEEE Computer Society[14] 上田和紀，細部博史，石井大輔 : ハイブリッド制約
Press, 1996, pp. 278–292.
言語 HydLa の宣言的意味論， コンピュータソフトウェ[ 7 ] Andrade, V. A. M., Comda, L. D. J. and Stolﬁ,
ア，Vol.28, No.1, 2011, pp. 306–311.
J.: Aﬃne Arithmetic, INTERVAL’ 94, St. peters[15] Wolfram Research, Inc., Mathematica,burg, Russia, March, 1994.
http://www.wolfram.co.jp/products/
http://www.dcc.unicamp.br/˜stolﬁ/EXPORT/bibliography/
mathematica/index.html.