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

パラメタを含むハイブリッドシステムに対するアフィン演算を用いた記号シミュレーション
松本 翔太 上田 和紀
本研究の目的は，離散変化と連続変化の両方をともなう動的システムであるハイブリッドシステムの高信頼シミュ
レータの構築である．特に，パラメタを含むハイブリッドシステムの厳密なシミュレーションを行うことで，単なるモデルの到達可能範囲の解析だけではなく，望ましい性質を満たすためのパラメータ設計に役立てることを目標とす
る．区間演算を用いたシミュレータは到達可能範囲の包囲を求めることができるが，モデル中の不確定量がモデルの振る舞いに与える影響を計算することは依然として挑戦的課題である．
本論文では，記号計算を区間ニュートン法およびアフィン演算と連携させたシミュレーションアルゴリズムを提案する．本手法では，不確定量間の 1 次の依存関係を保存したまま計算を進めることが可能である．さらに，著者らが開発してきた記号シミュレータに本アルゴリズムを実装し，その性能を例題を用いて評価する．

The purpose of this research is to develop a highly reliable simulator of hybrid systems, i.e., systems involving
both discrete change and continuous evolution. In particular, we aim at rigorous simulation of parametric hybrid systems, which enables not only the analysis of model’s possible behavior but also the designof parameters that realize desired properties. Simulators with interval arithmetic can reliably compute areachable set of states, but preserving the dependency of uncertain quantities in models is still challenging.
In this paper, we discuss a simulation method that is based on symbolic computation and cooperates
with the interval Newton method and affine arithmetic, which is able to preserve first-order dependency ofuncertain quantities. We implemented the algorithm on the symbolic simulator we have been developing
and evaluated the performance of the method with example models.

とんどのツールは区間演算 [15] に基づいており，モデ1 はじめにルの到達可能状態の過大近似（over-approximation）

ハイブリッドシステムは離散変化と連続変化の両方を求めることができる．

をともなう動的システムである．ハイブリッドシステ

感度解析やコントローラの設計などにおいて，モデ

ムのシミュレーションにおいては，計算機上の誤差がルに含まれる不確定パラメータがそのふるまいに与

定性的に誤った結果を生むことにつながるため，モデえる影響の計算はさまざまな応用が考えられる．しかルの設計の正しさを保証するためには正しさが厳密にし，既存のツールではそれらの依存関係は無視される保証された計算技術を用いることが重要になってくる．

か，計算の途中で暗黙的に利用されるだけで，最終結

また，モデル化誤差や経年劣化など，モデルに含ま果には含められていなかった．本研究の達成目標は，れる不確定性を扱えることも重要である．ハイブリッ

パラメタを記号的に扱い，その依存関係を明示的に扱

ドシステムの厳密な解析を行うツールはこれまでにうことができるシミュレータを開発することである．

多く開発されてきている [4][6] [8] [17] [20]．その中のほ

この目標を達成するための課題として，(i) パラメタ

Symbolic Simulation of Parametric Hybrid Systems
with Affine Arithmetic
Shota Matsumoto, 早稲田大学基幹理工学部情報理工学
専攻, Waseda University.
Kazunori Ueda, 早稲田大学理工学術院情報理工学科,
Waseda University.

を含む微分方程式を解くことと (ii) パラメタを含むハ
イブリッドシステムの離散変化を厳密に計算することの 2 点があげられるが，本稿では後者に主眼を置き，
アフィン演算 [2] を用いた記号シミュレーションアル
ゴリズムについて論じる. 前者に関しては，微分方程式が線形であるか，もしくはパラメタを含む線形な微水槽の水がパイプを通して 2 つ目の水槽に流れ込ん

分方程式によって近似可能であることを仮定している

でいる．x1 と x2 は各水槽の水位を示しており，v1（近似誤差の包含に，新たなパラメタを用いることでと v2 は各水槽のバルブの状態を示している．vi = 1

対応する）．アフィン演算は区間演算を拡張した演算(0) はそれぞれバルブが開いている（閉じている）こであり，記号パラメータを用いることで不確定値間のとを示している．HydLa における変数はすべて（暗一次の依存関係を保存可能な手法である．

黙的に）時刻に関する関数変数である．本プログラム

ハイブリッドシステムのシミュレーションにおいては，時刻 0 において成り立つ性質を記述したいくつは，各離散変化が発生する時刻を求める必要があるかの制約から成り立っている．制約 INIT はシステム

が，その条件は解析解を持たない非線形方程式で記の初期状態を記述している．本モデルでは，1 つ目の述されている場合がある．このため，離散変化時刻の

水槽の初期水位は 1.9 ≤ x1 ≤ 1.9001 を満たす不定計算に関しても区間演算の技術を用いる必要がある．

値となっている. 制約 X1 と制約 X2 は各水位の連続

著者らの経験では，システムの連続挙動を解析解を持的挙動に関する制約である．時相演算子 [] はその制

つ線形微分方程式に限定しても，離散変化条件を表す

約が有効になった時刻以降常に [] のついた制約が成方程式がパラメタを含むために，しばしば記号的に立することを意味する．x1’ は時間微分 ∂x1/∂t を表解くことに失敗するケースに遭遇している．これをす. 含意記号 => の右辺は，その左辺（ガードと呼ぶ）

解決するため，本稿では区間ニュートン法と中間値のが満たされる時に有効になる．制約 V1 CONST と制定理を利用した離散変化時刻の計算方法を提案する．

約 V2 CONST は v1 と v2 が変化しないことを記述しこれらの技術を連携させることで，アフィン演算単体ており，これがこのモデルの通常時のふるまいとなっでは離散変化の計算に失敗するようなモデルに対している．制約 V1 OFF2ON，V1 ON2OFF，V1V2 OFF2ON，てもシミュレーションを行うことが可能になった．本V2 ON2OFF はバルブの開閉に関する制約であり，x1-

研究では提案手法を，著者らが開発してきたハイブのように変数のあとにマイナス記号がついたものはリッド制約言語 HydLa の記号シミュレータ上に実装
し，例題を用いた性能評価を行った †1 ．

左極限値 limtl ↑t x1(tl ) を表す．最後の 4 行はこれらの制約間の優先順位を << を用いて宣言しており，水位が閾値をまたぐ時以外はバルブの状態は変化しな

2 ハイブリッド制約言語 HydLa

いことを表現している．なお，より優先度の高い制約

本稿で提案するシミュレーションアルゴリズムへを持たない制約（本プログラムでは Vi CONST 以外のの入力は HydLa [19] プログラムとして与えられる．

制約）は常に有効な制約である．HydLa プログラム

HydLa は制約に基づいたハイブリッドシステムのモの宣言的意味はプログラムに記述された仕様を満た

デリングを行う．制約による記述は宣言的でありながす軌道の集合であり，各時刻における仕様は制約間のら，同期や条件分岐を含む制御構造も表現可能であ優先順位に違反しないもののうち，極大無矛盾な制約

る．さらに，制約を用いることで不確定性を自然にの部分集合として与えられる．

扱うことも可能であり，記号パラメータを用いたプログラムの記号実行に適している．本節では，例題を3 区間演算とアフィン演算

用いて HydLa の概要を紹介する．HydLa の構文や

区間演算 [15] は計算結果の取りうる範囲を保証す

意味論に関する詳細は [19] を参照されたい．図 1 はるための演算である．区間演算における各計算は実数

例題モデル（ [5] より引用）の概要を示しており，図の閉区間に対して定義されており， [l, u] という記法

2 は本モデルに対応する HydLa プログラムを示してで実数の集合 {x ∈ R | l ≤ x ≤ u} を表す. 本稿では，いる．本モデルは 2 つの水槽を含んでおり，1 つ目の実数区間のなす全体集合を I と記述する. 区間演算で†1 本稿は， [13] で発表予定の論文をベースにしたものである

は，計算誤差や不確定量の範囲をともに区間として表現することができる．通常の区間演算ではラッピン

れる．

x = x0 + x1 ǫ1 + · · · + xn ǫn
ここで，ǫi はノイズ変数と呼ばれる記号パラメタ
であり，[−1, 1] の範囲の値を取りうる不確定値である．x0 は x の中心値を表しており，xi (i > 0) は
各ノイズ変数の x に対する影響の大きさを表している．ǫi が 2 つ以上の量の中に現れた場合，それらの量の間には依存関係が存在する．アフィン形式は，
P
Pn[x0 − ni=1 xi , x0 +i=1 xi ] を計算することで対応する区間に直すことができる. 計算機上でアフィン演
図1

2 つの水槽のモデル算を行う場合，加算や乗算などの各計算を行うたびに，丸め誤差や非線形演算の近似誤差を包含するための新たなノイズ変数が導入される．ノイズ変数はINIT <=>
1.9 <= x1 <= 1.9001 /\ x2 = 1/\ v1 = 0 /\ v2 = 1.
X1 <=>[]((v1 = 0 => x1’ = -x1 - 2)/\ (v1 = 1 => x1’ = -x1 + 3)).
X2 <=>[]((v2 = 0 => x2’ = x1)/\ (v2 = 1 => x2’ = x1 - x2 - 5)).
V1_CONST <=> [](v1’ = 0).
V2_CONST <=> [](v2’ = 0).
V1_OFF2ON <=>[](v1- = 0 /\ x1- = -1 => v1 = 1).
V1_ON2OFF <=>[](v1- = 1 /\ v2- = 1 /\ x1- = 1 => v1 = 0).
V1V2_OFF2ON <=>[](v2- = 0 /\ x2- = 1 => v2 = 1 /\ v1 = 0).
V2_ON2OFF <=>[](v2- = 1 /\ x2- = 0 => v2 = 0).

計算結果の精度を向上させるために有用ではあるが，使用されるノイズ変数の数が増えるにつれて計算コストも大きくなってしまう．そうした精度とコストのトレードオフを制御するために，ノイズ変数の数を削減するアルゴリズムが提案されている [9]．この削減
アルゴリズムはノイズ変数のうち，削減した場合の精度の低下が小さい順に指定された数までノイズ変数を減らすことが可能である．本稿で用いる手法は，各離散変化の計算の際にこの削減アルゴリズムを利用している（2 節）．

4 シミュレーションアルゴリズム
本節で提案するアルゴリズムは与えられた HydLa
プログラムの軌道を記号的に計算するものである．図

INIT, X1, X2,(V1_CONST, V2_CONST)<< (V1_OFF2ON, V1_ON2OFF,
V1V2_OFF2ON, V2_ON2OFF).
図2

2 つの水槽を記述した HydLa プログラム

3 にアルゴリズムの概要を示す．このアルゴリズム
は基本 HydLa [19] のプログラムを入力とし，軌道の集合を出力する．基本 HydLa のプログラムは，元のHydLa プログラム中の個々の制約間の優先度に関する仕様を，その仕様を満たす制約の部分集合を要素とグ効果や数値間の依存関係の無視により，計算結果のする半順序集合に変換することで得られる（この半順区間の幅が広がってしまうという欠点がある．この欠序集合中の順序関係として，部分集合間の包含関係を点を解決するための手法は数多く研究されており [3]用いる）．別の言い方をすると，制約集合 S が制約[15]，その 1 つとしてアフィン演算 [2] が存在する．ア

C を含むならば，S は必ず C より高い優先順位を持フィン演算は区間の 1 次の依存関係に着目した区間

つような制約集合のみを，半順序集合中の要素として演算の拡張である．アフィン演算では，それぞれの数

選ぶ．

量 x はアフィン形式と呼ばれる以下の形式で与えら

本アルゴリズムは 2 つのフェーズを交互に切り替えながら進行する．1 つ目のフェーズは離散変化を扱う

Point Phase であり，2 つ目のフェーズは連続変化を扱う Interval Phase である．各フェーズでは，HydLaの意味論に従い制約の極大無矛盾集合（MCS ）を計算して S に代入し，記号パラメータに関する条件をP に代入する (10 行目–16 行目)．MCS は各フェーズに対応した無矛盾性判定を行う関数を引数とする
高階関数である．9 行目の Subst は S 中の各数式に現れる時刻変数 t に対し，現在時刻 T を代入する．

Point Phase の最後に，その時刻で値が不確定な変数に対して記号パラメータを導入し，その値の包囲を

Enclose 関数によって計算する（14 行目）. Enclose
関数はアフィン演算を用いており，最終的に残すノ
イズ変数の数を指定することができる．6 節において，ノイズ変数の数が性能に与える影響について比較を行う．本アルゴリズムでは，アフィン演算との親和性を高めるため，P 中の全ての記号パラメタ（piとする）は pi ∈ [−1, 1] を満たすように正規化されている．例えば，下限と上限がそれぞれ x と x で表される不確定値は，(x + x)/2 + (x − x)/2 × pi のようにする．Interval Phase では，再び MCS を求めるとともに，成立したガードの集合 A+ と成立しなかったガードの集合 A− を得る（16 行目）. その後，

S 中の微分方程式を解き（17 行目）次の離散変化時刻の計算を行う．
（21 行目）．Subst は微分方程式の解を S と g と ¬g 中の各数式に対して代入している．一般に，離散変化時刻の候補はパラメタの値の不確定性によりが複数存在しうる．例えば，x の軌

入力: HydLa: 基本 HydLa のプログラム,
MaxT : 最大シミュレーション時間
1: MS := TopologicalSort(HydLa)// 候補となる制約集合のリスト
2: V := GetVariables(HydLa)
3: T := 0// 現在時刻
4: S := true // 現在成立している制約集合
5: P := true // 記号パラメタの条件
6: Gp := ∅ // 直前の離散変化の原因となったガードの集合
7: while T < MaxT do
8:// Point Phase (PP)9:
S := Subst(S, T )10:(S, P, Gp , , ) :=MCS (S, MS, P, T, Gp , CheckConsistencyPP)
11:
if S = false then
12:
break
13:
end if
14:(S, P) := Enclose(AddParameters(S, P, V ))15:// Interval Phase (IP)16:
(S, P, Gp , A− , A+ ) :=MCS (S, MS, P, T, Gp , CheckConsistencyIP)
17:
S := SolveDifferentialEquation(S)18:
if S = false then
19:
break
20:
end if
21:(MinT,P,Gp ) := CompareMinTime({FindMinTime(Subst(g, S), P, Gp ) | (g ⇒ c) ∈A− }∪ {FindMinTime(Subst(¬g, S), P, Gp )) | (g ⇒c) ∈ A+ }∪ {(MaxT −T , true)})22:
T := MinT + T
23: end while
図3

シミュレーションアルゴリズムの概要

道が x(t) = −(t − 1)2 + x(0)，−1.5 ≤ x(0) ≤ 1.5
であり，離散変化条件が x(t) = 1 ∨ x(t) = −10 で場合に分岐させる必要がある．著者らがこれまでに開ある場合には，2 つの候補時刻が存在する．第 1 のp
候補は 1 ≤ x(0) ≤ 1.5 の場合の t = 1 − x(0) − 1

発した記号シミュレータ [14] の記号計算はそうした分であり，第 2 の候補は −1.5 ≤ x(0) < 1 の場合のpt = 1 − x(0) − 9 である. このように複数の候補がている区間ニュートン法はそのような場合を扱うこと現れるのは，パラメタの条件がある種のコーナーケー

アルゴリズムに拡張することを考えている．なお，本ス（例えばボールが壁にちょうど接するような軌道を稿で扱う例題モデルはこうした分岐が発生しないモ描く場合）を含む場合である．CompareMinTime やデルである．本アルゴリズム中の FindMinTime とFindMinTime では，複数の候補が存在するかどうか

Enclose 以外の手続きはすべて記号的に計算され，すを判定することが可能である．もし複数の候補が存在べての不確定量は P 中のパラメタとして扱われる．

する場合，取りうるすべてのふるまいを計算するためにはパラメタ空間を分割し，以降の計算をそれぞれの岐に対応することができるが，FindMinTime で用いができないため，今後の課題として分岐に対応した

4. 1 FindMinTime
FindMinTime 関数は離散変化時刻を計算する関数
であり，離散変化時刻としてパラメタを含むことを許している．HydLa プログラムにおいて，離散変化はプログラム中のガードの成否が変化した時に発生す
る．ガードは連立不等式および方程式によって記述され，時刻に関する制約条件とみなすことができる．このため，FindMinTime の目的は与えられた条件の成否が変化する（成立から不成立へ，もしくはその逆）
最も早い時刻を計算することとなる．
図 4 に FindMinTime のアルゴリズムを示す. まず，与えられた条件 G から原子境界条件のリストを得る.
原子境界条件は，G 中の各原子条件（方程式もしくは不等式）の関係演算子を等号に置き換えることで得られる．例えば，GetAtomicBoundaryConditions(t >0 ∧ t2 = 1 ∧ t ≤ 1) の結果は {t = 0, t2 = 1, t = 1} となる．次に，SolveTimeEquation によって各方程式の解区間を計算する．SolveTimeEquation は与えられた方程式を解析的に解くことができるなら，その解

入力: G: 微分方程式の解をガードに代入してえられた，時刻 t に関する条件，
P : パラメタの条件
Gp : 直前の離散変化の原因となったガードのリスト
出力: tmin : ガードの成否が変化する最も早い時刻,
glist : tmin において成否が変化する原子ガードのリスト
1: glist := GetAtomicBoundaryConditions(G)2: tlist := ∅; Mapg := ∅3: for gb ∈ glist do
4:(sols, P ) := SolveTimeEquation(gb , P, Gp )5:
for sol ∈ sols do
6:tlist .add(sol, gb )7:
end for8:
Mapg .insert(gb , ConsistentAtInitialTime(gb , Gp ))9: end for
10: while tlist 6= ∅ do
11:
(tmin , glist ) := PopMinimumTime(tlist , P )12:(Mapg , satisfied) :=CheckAndUpdateGuards(Mapg , glist , G)13:
if satisfied = true then
14:
break
15:
end if
16: end while

析解を返す．そうでなければ，方程式を 4. 2 節に示

図4

FindMinTime のアルゴリズム

す方法で解き，解のリスト sols と，解区間に対応する新たなパラメータの条件を含んだ P を返す．それ定を行う．g が境界を通りすぎているため，更新時のぞれの解は，g と組にして tlist に入れておく．これ

判定が 3 行目から 9 行目のものと異なることに注意にともない，各ガードが現在の Interva Phase の初期

されたい．

時刻で成り立っているかどうかの対応表として Mapgを作成しておく（8 行目）. 9 行目から 16 行目では，

4. 2 SolveTimeEquation

tlist 中の各時刻を始点とした時刻区間について，ガー

SolveTimeEquation は解析解を計算できない場合，

ド全体の成否を調べていく． PopMinimumTime は

精度保証数値計算による計算を行う．この手続きは 2

tlist から最も早い時刻 tmin を取り出し，tmin におい

つの段階に分けられる．

て成否が変化する原子境界条件 glist とともに返す．

まず，与えられた方程式 gb を区間ニュートン法

CheckAndUpdateGuards は tmin から始まる時刻区[15] によって解く. 区間ニュートン法では，初期区間間における，ガード全体の成否と，各ガードについて

X (0) から始めて，
X (k+1) = X (k) ∩ N (X (k) )成否を更新した対応表 Mapg を返す.
図 5 に CheckAndUpdateGuards のアルゴリズムを示す．3 行目から 9 行目では， tmin ちょうどにおという漸化式によって区間を狭めていく．なお，N (X)はける各原子条件の成否を計算している．関係演算子

N (X) = m(X) − f (m(X))/f ′ (X)（relop(g)）が {‘=’, ‘≤’, ‘≥’} のいずれかの場合は成であり，m(X) は X の中点である．X (0) が真の解を立し，そうでなければ成立しない．10 行目ではガー含んでいれば X (k) の区間幅は 2 次収束することが証

ド全体の成否を調べている． 13 行目から 24 行目で明されている．この計算は X (0) が gb の複数の解をは，tmin から始まる時刻の開区間において同様の判含んでいる場合には分岐し，それぞれが真の解を含む

入力: Mapg : 直前の時刻区間における，ガードとその成否の対応表，
glist : 現在時刻で成否が変化する原子ガードのリスト，
G: 全体のガード
出力: Mapg : 現在の時刻区間における，ガードと成否の対
応表，
satisfied: 現在の時刻区間における G の成否
1: Mapprev := Mapg
2: satisfied := false
3: for g ∈ glist do
4:
if relop(g) ∈ {‘=’, ‘≤’, ‘≥’} then5:
Mapg .replace(g, true)
6:
else
7:
Mapg .replace(g, false)8:
end if
9: end for
10: if G.satisfiedBy(Mapg ) then11:satisfied := true
12: end if
13: for g ∈ glist do
14:
if relop(g) = ‘=’ then15:
Mapg .replace.(g, false)
16:
else if relop(g) = ‘6=’ then17:
Mapg .replace.(g, true)
18:
else
19:
Mapg .replace.(g, ¬Mapprev .getV alue(g))
20:
end if
21: end for
22: if G.satisfiedBy(Mapg ) then23:
satisfied := true
24: end if

メタの 1 次の依存関係を保存した記号解を計算することを目的とする．このステップの仕様は次の通りである：
入力: f (t, p~) : R × [−1, 1]n → R,

T: 区間ニュートン法によって計算された時刻区間
出力: Tresult (~p) : [−1, 1]n → I （方程式 f (t, p~) = 0の解 t(~p) : [−1, 1]n → R を包含する，パラメタ化
された時刻の式）
以下，Tresult (~p) の計算手続きについて述べる．ま~ ∈ In ，t ∈ T ，p
~ に対し，中間値の定理より
ずP~∈P
f (t, ~p) ∈ f (Tm , P~m )~ ) ∂f (T, P~)∂f (T, P∂f (T, P~ )
+(
,...,
)
,
∂t
∂p1
∂pn
~m1 , . . . , pn − P~mn )· (t − Tm , p1 − P(1)~m は各要素がが成り立つ．Tm は T の中点であり，PP~ の対応する要素の中点であるようなベクトルである．以下は式（1）の証明である．

Proof. 中間値の定理より，連続かつ微分可能な関数
h : Rn → R と閉区間 [a, b] ∈ In に対し，以下を満たす c ∈ [a, b] が存在する．

h(b) = h(a) + ∇h(c) · (b − a)
ここで，[a, b] ⊆ I を満たすような区間 I を考えると，

I は h(c) ∈ h([a, b]) ⊆ h(I) を満たす．これにより，以下が言える．

図5

CheckAndUpdateGuards のアルゴリズム

h(b) ∈ h(a) + ∇h(I) · (b − a)
上 式 に お い て ，h，b，a，I を そ れ ぞ れ f ，(t, p~)，ような複数の区間を計算することができる．
なお，もし gb が Gp に含まれている場合，gb がは
現在の時刻区間の始点において成立するため，区間ニュートン法は始点を含む時刻区間に対しても分岐してしまう．そのような分岐は無視する必要があり，この判定を行うために SolveTimeEquation は引数として Gp を取っている．
次に，区間ニュートン法で得られた結果を用いてパラメタの 1 次項を保存した解を計算する．区間ニュートン法の結果は，gb 中のパラメタを区間に置き換えて計算することで得られた単なる区間解であり，パラメタの依存関係は放棄されていることに注意されたい．このステップでは，中間値の定理に基づいてパラ~m )，(T, P~ ) に置き換えることで，式（1）を(Tm , P
得る.
ここで各パラメタの範囲は [−1, 1] に正規化してあ~m は ~0 と見なすことができ式（1）は次のるため，Pように簡約化できる．
f (t, p~) ∈f (Tm , ~0)~)~ ) ∂f (T, P~ )∂f (T, P∂f (T, P
,...,
)
,
+(
∂t
∂p1
∂pn
· (t − Tm , p1 , . . . , pn )(2)式（2）の右辺を fmean とし，方程式 fmean = 0 を tについて解くことで次式を得る．
(f∂p1 , . . . , f∂pn )f (Tm , ~0)t=−·p~ + Tm −
f∂t
f∂t

(3)

~ )/∂t とこ こ で ，f∂t と f∂pi は そ れ ぞ れ ∂f (T, P~ )/∂pi を表す区間値である. これらの区間∂f (T, P

5 提案手法の実装

を記号パラメタとして扱うと，計算コストは急速に増

2 節に示した手法を，著者らが開発を続けている

大する可能性がある．これを避けるため，各離散変化

記号シミュレータ HyLaGI [14] [18] 上に実装した．

時刻の計算において新たに導入する記号パラメタを 1

HyLaGI は C++で実装されており，現在は Mathe-~ を満たすすべてのつに限定することを考える．~
p∈P

matica をバックエンドの制約ソルバとして利用してp~ について，次の性質が成り立つ．
f∂p1 , . . . , f∂pn
(f∂p1 , . . . , f∂pn )·p~ ∈ mid(
)·p
~
f∂t
f∂t
n
X
f∂pi
)
+ [−1, 1] ×
rad(
f∂t
i=1
(4)

いる．シミュレーションアルゴリズム全体の流れは式（3）と式（4）より，元の方程式 f (t, p~) = 0 を包性能は理想的とは言い難いが，記号的手法と精度保証

含する記号解として次の Tresult (~p) が得られる．

数値計算技術を組み合わせることが容易な，柔軟なプf∂p1 , . . . , f∂pn
Tresult (~p) = mid(
)·p
~
f∂t
n
X
f∂pi
f (Tm , ~0))−+ Tm + [−1, 1] ×
rad(
f∂t
f∂t
i=1
(5)
Pn
~
ここで，[−1, 1] ×rad(f∂pi /f∂t ) − f (Tm , 0)/f∂t
i=1

は 1 つの区間に還元することができ，そうすると式（5）の右辺はアフィン形式と見なせる．この還元後のTresult (~p) がこの手続きの最終結果となる．
Example 1. SolveTimeEquation に対し，以下の入力を与える．
gb := 2 − exp(t + py /1000) + px /1000 = 0,

P := −1 ≤ px ≤ 1 ∧ −1 ≤ py ≤ 1,
Gp := ∅
まず，区間ニュートン法で実数区間 [0.69163, 0.69466]を得る. 次に，Tresult (~p) を以下の計算によって得る．
Tresult ((px , py ))
f∂py
f∂px
= mid()px + mid(
)py + Tm
f∂t
f∂t
f∂py
f∂px
+ [−1, 1] × (rad() + rad(
))
f∂t
f∂t
~
f (Tm , 0)−f∂t
= 5.0000 × 10−4 × px
−4

− 9.9999 × 10

C++で実装されており，Mathematica は制約の連言の無矛盾性判定，常微分方程式の求解，離散変化時刻
に関する最小化問題，算術式や論理式の変形に用いている．6 節に示す通り，バックエンドの制約ソルバのラットフォームとして利用可能である．現在の実装の健全性はブラックボックスとして使用しているバック
エンドソルバの健全性にも依存している部分がある
が，今後は計算手続きをさらに細かくモジュール化していくことでブラックボックスとして使用する機能の粒度を小さくすることを課題としている．本実装において，区間演算やアフィン演算およびノイズ変数の削減に関しては KV ライブラリ [10] を利用している．

HyLaGI の出力は与えられた HydLa プログラムの仕様を満たす軌道の集合である．出力される軌道は各時刻の不確定な変数値に対し，記号パラメタを用いることで表現されている．例えば図 2 の例題プログ
ラムに対しては，x1 の初期値に対応する記号パラメ
タが導入される．区間演算における区間値や，アフィン演算におけるノイズ変数も，同様に記号パラメタとして扱われる．

6 実験結果
本節では，Enclose と FindMinTime において，異なる方法で計算を行った場合の結果を比較する．比較のために 2 つの例題を使用する．

× py
6. 1 2 つの水槽の例題+ 0.69314−5

1 つ目の例題は 2 節に示したものである. 取りうる

+ 1.0321 × 10 × [−1, 1].
なお，この例では話を簡単にするために仮数部を 10

軌道の 1 つを図 6 に示す．横軸は時刻であり，縦軸

進 5 桁に制限しているが，後述の実装は IEEE754 のがそれぞれの水位を示す．本モデルを，以下に示す 2

倍精度浮動小数演算に基づいたものになっている．

つの方法でシミュレーションした．

図7
図6

水位の区間幅

水位の軌道 （実線: x1, 破線: x2）

1. 区 間 ニュー ト ン 法 の 結 果 を そ の ま ま SolveTimeEquation の結果として用い，Enclosure の中でアフィン演算の代わりに区間演算を用いた（Newton&IA と略記）．

2. 節に示した通り，SolveTimeEquation で区間
ニュートン法と中間値の定理を利用し，Enclose
でアフィン演算を用いた（Mean&AA と略記）．なお，ノイズ変数の最大数は 6 個に制限する
図 7 は各フェーズにおける各変数値の区間幅の合計

図8

2 つの水槽の例題の実行時間を示しており，図 8 は各ステップの実行時間を示している．横軸はシミュレーションのステップ数を表し

6. 2 Bouncing Ball on a Sine Curveており，Point Phase と Interval Phase を合わせて 1

2 つ目の例題は，2 次元空間中において正弦波上のステップとしている．

床の上で質点が跳ねるモデルである．図 9 にプログ

図 7 において，最初の区間幅は 0.0001 である．こラムを示す．本プログラムでは，x と y が質点の位置のプログラムでは，水位の理想的なふるまいは初期水を示し，e が反発係数を示す．cont，s，c は質点の位に不確定性があっても一定のふるまいに収束してバウンドを記述するための補助変数である．図 10 にいく．Newton&IA では不確定値の依存関係を無視し本モデルの軌道の例を示す．

ているため，区間幅が収束せず，61 ステップ後にシ

図 11 に各フェーズの区間幅を，図 12 に実行時間

ミュレーションが破綻する．一方，Mean&AA では依存を示す．本モデルは初期状態に不確定性を持たず，関係を扱うことに成功しており，区間幅が 10−7 未満

区間幅の増大は純粋に計算誤差によって発生する．

に収束する.

Newton&IA と Mean&AA の意味は 6. 1 節のものと等し

図 8 に示す通り，Newton&IA の計算時間は Mean&AA

い. 本実験ではそれに加え，保存するノイズ変数の数より短い．どちらの方法でも計算時間はステップ数にを 5，9，13 と変化させる（それぞれ d5，d9，d13 と対して定数オーダーに収まっており，これは Enclose略記する）．本モデルにおいては，x, x’, y, y’，t のによって数式の精度保証近似を行っているためと考え

5 つの変数に対してアフィン形式を用いるため，ノイられる．もし Enclose を用いない場合，変数値を表

ズ変数の最小数が 5 であることに注意されたい．

す数式のサイズは増大し続け，実行時間はこのプログ
ラムの場合線形に増大していくことになる．

図 11 を見ると，どの条件でシミュレーションを行っても区間幅はステップ数に対して指数的に増大することがわかる．また，ノイズ変数を増やすことで増大速度が減少することもわかる．区間幅が増大して正弦波

INIT <=>
x = 0 /\ x’ = 0 /\ y = 10 /\ y’ = 0/\ e = 1 /\ [](e’ = 0).
FALL <=> [](cont = 1 => y’’ = -10).
CONSTX <=> [](cont = 1 => x’’ = 0).
SC <=>[](s = cos(x-)/(1 + cos(x-)^2)^(1/2)/\ c = 1 /(1 + cos(x-)^2)^(1/2)).
BOUNCE <=>[]( y- = sin(x-) => cont = 0/\ x’ = ((-e) * s-^2 + c-^2) * x’+ ((e+1) * s- * c-) * y’/\ y’ = ((e+1) * s- * c-) * x’+ (s-^2 + (-e) * c-^2) * y’-).

図 11

跳ねる質点の区間幅

図 12

跳ねる質点の実行時間

INIT, SC, FALL, CONSTX,[](cont = 1) << BOUNCE.
図9

図 10

跳ねる質点に対応する HydLa プログラム

跳ねる質点の軌道（実線：床，破線：質点）の範囲と時刻に関して branch and prune アルゴリズの振幅を超えると，FindMinTime は失敗し，それ以ムを適用する，domain contraction を利用している．

降のシミュレーションを行うことはできなくなる．

Flow* は 6 節に示した例題の区間包囲を高速に計算

図 12 を見ると，ノイズ変数を増やすと実行時間もすることができ，区間幅も本手法と同様に収束する

増えるが，各ステップの実行時間には何らかの上限が，ガード部やインヴァリアント部に非代数的関数があることがわかり，一定以上には増大していない．

を記述できないため，6. 2 節の例題は扱うことができ

また水槽の例題と同じく，Newton&IA の計算時間はない．

Mean&AA よりも少ない．

Acumen [20] は独自の手続き型言語を入力とするハ
イブリッドシステムのシミュレータである．3D オブ

7 Related Work

ジェクトの描画を含む豊富なグラフィック機能に加え，

本節では，ハイブリッドシステムのふるまいを厳密

Zeno 時刻以降の軌道の包囲も計算できる精度保証シに計算するツールを関連研究として紹介する． Flow*

ミュレーションをサポートしている [1]．著者らが実[4] はハイブリッドオートマトン [7] の到達可能範囲を験したところ，Acumen の現在の実装は 2 つの水槽計算するツールである． Flow*では Taylor model にの例題を扱うことができたが，計算結果の区間幅は収基づき，非線形微分方程式の解包囲を計算すること

束せず，12 回バルブを開閉したところでシミュレーができる．離散変化の計算に関しては，変数の初期値

ションが進まなくなった．

dReach [11] は常微分方程式を扱うことができる

る．これは，区間ニュートン法における N (X) の計

背景理論付き SAT（Satisfiability Modulo Theories,

算において，f (m(X)) と f ′ (X) の両方が 0 を含むとSMT) ソルバ dReal をベースとした，ハイブリッド

計算に失敗するためである．この問題を解決するたシステムの有界モデル検査器である．dReach はハイめ，パラメタ空間を分割して計算を進めることを計画

ブリッドオートマトンの有限ステップの実行を展開している．また，解析解を持たない非線形微分方程式

し，SMT の問題に帰着する．dReach は区間制約伝を，パラメタを含む線形微分方程式によって包囲する

搬に基づいており，与えられた性質の充足可能性を判ことや，関連ツールとの詳細な性能比較も重要な課題

定することができる．dReach は充足可能性の証拠とである．

して軌道の包囲も出力することができるが，その結果はパラメタを含んでいない．

KeYmaera [17] および KeYmaera X [16] は differential dynamical logic に対する定理証明器であり，謝辞 精度保証数値計算技術に関して，貴重なご意見をいただいた柏木雅英教授に感謝いたします．本研究の一部は, 科学研究費補助金 (26280024, 15K12010)の補助を得て行った.

記号計算を用いてハイブリッドシステムを扱うツールという点で本研究と関連している．KeYmaera は

参 考 文 献

定理証明器としてのアプローチを取っているため，無[ 1 ] Duracz, A., Bartha, F.A., Taha, W.: Accurate
rigorous simulation should be possible for good designs, International Workshop on Symbolic and Numerical Methods for Reachability Analysis, 2016,
pp. 1–10.
[ 2 ] de Figueiredo, L. H. and Stolfi, J. : Affine Arithmetic: Concepts and Applications, Numerical Algorithms, Vol. 37, Issue 1-4, 2004, pp. 147–158.
[ 3 ] Berz, M. and Makino, K. : Higher order multivariate automatic differentiation and validated computation of remainder bounds, WSEAS Transactions on Mathematics, 1998, Vol. 3, Issue 1, pp. 37–
44.
[ 4 ] Chen, X., Abraham, E. and Sankaranarayanan,
S. : Flow* : An Analyzer for Non-Linear Hybrid
Systems. CAV 2013, LNCS 8044, Springer-Verlag,
2013, pp. 258–263.
[ 5 ] Chen, X., Schupp, S., Makhlouf, I. B., Abraham,
E., Frehse, G., Kowalewski, S. : A Benchmark Suite
for Hybrid Systems Reachability Analysis, In 7th
NASA Formal Methods Symposium, Springer, 2015
pp. 408–414.
[ 6 ] Frehse, G., Guernic, C. L., Donzé, A., Cotton,
S., Ray, R., Lebeltel, Olivier., Ripado, R., Girard,A., Dang, T., Maler, O. : SpaceEx : Scalable Verification of Hybrid Systems, CAV 2011, LNCS 6806,
Springer-Verlag, 2011, pp. 379–395.
[ 7 ] Henzinger, T. : The Theory of Hybrid Automata, LICS’96, IEEE Computer Society Press,
1996, pp. 278–292.
[ 8 ] Hickey, T. J. and Wittenberg, D. K. : Rigorous Modeling of Hybrid Systems Using Interval
Arithmetic Constraints, HSCC 2004, LNCS 2993,
Springer-Verlag, 2004, pp. 402–416.
[ 9 ] Kashiwagi, M. : An algorithm to reduce thenumber of dummy variables in affine arithmetic,
15th GAMM-IMACS International Symposium on
Scientific Computing Computer Arithmetic and

限時間に関する性質も証明可能という長所があるが，証明手続きは本質的にユーザ対話型のものとなり，自動的に軌道を求めることができるシミュレータとは対照的な立場にある．KeYmaera はユーザのモデルに対する理解を支援するためにハイブリッドシステムのシミュレーション機能も有してはいるが，あくまで定理証明の補助としての機能であり，モデルに含まれる不確定性を扱ってはいない．また，KeYmaera や

KeYmaera X は完全に記号計算に基づいているという点でも HyLaGI と異なっている．

8 まとめと今後の課題
本稿ではパラメタを含むハイブリッドシステムに対するシミュレーションアルゴリズムを提案した．本ア
ルゴリズムは記号計算に対し，区間ニュートン法，アフィン演算，中間値の定理を用いた方程式求解という3 つの技術を連携させたものである．本手法を実装したシミュレータによりパラメタ間の 1 次の依存関係を含んだ軌道の式を計算することができ，計算結果精度
は通常の区間演算を用いたものより向上が見られた．
今後の課題として，広い範囲を持つパラメタの設計問題を解くために，不確定値が広い範囲を持つようなモデルを扱うことがあげられる．現在の Find-

MinTime は不確定値の範囲が狭く，パラメタの条件に関する場合分けが発生しないことを前提としていVerified Numerical Computations, 2012.
[10] 柏木 雅英 : kv library, http://verifiedby.me/kv/.
[11] Kong, S., Gao, S., Chen, W., Clarke, E.:dReach: δ-reachability analysis for hybrid systems,
TACAS 2015, LNCS 9035, 2015, pp. 200–205.
[12] Lunze, J. : Handbook of Hybrid Systems Control : Theory, Tools, Applications, Cambridge University Press, 2009.
[13] Matsumoto, S., Ueda, K.: HyLaGI: Symbolic
Simulation of Parametrized Hybrid Systems withAffine Arithmetic, 23rd International Symposiumon Temporal Representation and Reasoning (TIME
2016), 2016, accepted.
[14] Matsumoto, S., Kono, F., Kobayashi, T., Ueda,
K.: HyLaGI: Symbolic Implementation of a Hybrid Constraint Language HydLa, Electronic Notes
in Theoretical Computer Science, Vol. 317, 2015,
pp. 109–115.
[15] Moore R. E., Kearfott R.B., Cloud M. J.: Introduction to Interval Analysis, Society for Industrial

and Applied Mathematics, 2009.
[16] Fulton, N., Mitsch, S., Quesel, J.D., Völp, M.
and Platzer, A.: KeYmaera X: An axiomatic tactical theorem prover for hybrid systems. CADE 2015,
LNCS 9195, 2015, pp. 527–538.
[17] Plätzer, A., Quesel, J.D. : KeYmaera : A Hybrid
Theorem Prover for Hybrid Systems. IJCAR 2008,
LNCS 5195, Springer-Verlag, 2008, pp. 171–178.
[18] HydLa, http://www.ueda.info.waseda.ac.jp/hydla/.
[19] Ueda, K., Matsumoto, S., Takeguchi, A.,
Hosobe, H. and Ishii, D. : HydLa : A High-Level
Language for Hybrid Systems. Second Workshop on
Logics for System Analysis(LfSA 2012), 2012, pp. 3–
17.
[20] Zeng, Y., Rose, C., Brauner, P., Taha, W.,
Masood, J., Philippsen, R.., O’Malley, M. and
Cartwright, R. : Modeling Basic Aspects of CyberPhysical Systems, Part II. 11th IEEE International Conference on Embedded Software and Systems(ICESS 2014), 2014, pp. 550–557.