若槻祐彰　松本翔太　伊藤剛史　和田努　上田和紀時間経過に伴い状態が速続変化したり、状態や方程式系自体が離散変化したりするシステムのことをハイプリッドシステムと言う.  HydLaは, ハイプリッドシステムを制約階層によってモデリングする宜言型言語である.  HyLaGIは数式処理に基づいたHydLaの処理系であり、記号パラメータを含むシステムを扱える。例えばボール同士の衝突における三体同時衝突のような、特異点の近傍の運動に興味がある場合に、特異点から微小値εだけずれた運動を計算し、εを限りなく0に近づけることで解析を行うことができる．HyLaGIはそのような解析をサポートする機能を有しており、特異点の解析以外にも誤差の解析に利用することも可能である。本研究では本機能を利用していくつかの例題の解析を行うことで，その有用性を示す. Hybrid systems are dynamical systems that are composed of continuous behaviors and discrete behaviors. HydLa is a constraint-based modeling language for Hybrid systems.  HyLaGI is an implementation of HydLa. HyLaGI can simulate systems that have simbolic parameters. For example, when we are interested in move- ments that are close to a singularly point, we can simulate movements that are away from the singular point by a symbolic parameter ε and then analyze movements by taking the limit of ε . HyLaGI supports such simulation using ε . In this paper, we show that this function is useful for analyzing a number of examples.
1研究背景
1.1ハイブリッドシステム
ハイプリッドシステム[3]とは，時間の経過に伴い状態が連続変化したり，状態や方程式系自体が離散変化したりするシステムのことである。物理現象やコンピュータ調御されたシステムの多くは，ハイブリッドシステムとしてみなすことができる。このように物理学だけでなく，調御工学や生命工学などのさまざまな分野においてシステムをモデリングする際に利用することができる.
1.2 ハイブリッド制約言語HydLa
HydLa[6][5]とはハイプリッドシステムをモデリングする宣言型言語である。システムの挙動を方程式や論理式を用いた制約によって記述する。制約間に優先順位をつけることで，システムを簡潔に表現することができろ.例えば，地面を跳ねる質点をモデリングすると図1のように装すことができる. 1 2 3
4図１　床で跳ねる質点のプログラムここで，変数ｙは質点の位置を表している． INITは，質点の初期位置が10，初速度が0であることをを表す制約である。FALLは，時刻0以降において質点の加速度が−10であることを表す制約である. BOUNCEは，時刻0以降において質点の位置の左極限値が0の場合，質点の速度は現在の速度の左極限値に−4/5をかけた値である，つまり質点が地面に衝突した時に，速度が減衰して跳ね返ることを表している制約である．この=>記号の左辺をガード条件といって，この論理式が真のときに右辺の方程式が有効になる制約を表している．ここで地面に衝突した瞬間は，FALLとBOUNCEの二つの制約が矛盾する．そのような場合にどちらの制約を優先するのかを制約階層[1]を用いて表したのが，4行目のINIT，FALL << BOUNCEである。ここではFALLとBOUNCEが矛盾した場合BOUNCEを優先するように記述してある.
1. 3    HydLa処理系HyLaGI
HyLaGI[4]は数式処理を用いたHydLaの処理系であり，浮動小数点誤差のないシミュレーションが可能である。また記号実行が可能で，システム内の不定値をパラメータとしてそのまま扱うことができる.例えば，図１をシミュレーションした結果として図2のように出力する．3,20行目のPPは離散変化するフェーズ(Point Phase)を表し，11行目のIPは連続変化するフェーズ(Interval Phase)を表している.それぞれのフェーズの出力内容は，その時に採用されなかった制約(unadopted modules)や衝突した制約(unsat mod, unsat cons),成立(positive)もしくは不成立(negative)に変化したガード条件，またその時の時刻(t)と各変数(y,y',y'')の値となっている.
2　誤差を用いたモデル解析
HydLaとHyLaGIを用いることで，物理現象やコンピュータ制御されたシステムのふるまいをシミュレーションすることができる。しかし，現実にモデルのふるまいを考えると与える初期値に誤差が入ることが考えられる。この場合，誤差の影響により，システムのふるまいが誤差を含まない場合のふるまいと大きくことなる場合がある。そのためシステムの性質を解析する際に，モデルに誤差が与える影響を調べることは大変亜要である．このような誤差解析だけでなく，ある値とその近傍を対象としたシミュレーションや無限大や無限小を扱う場合にも微小な誤差は利用可能である。本研究では微小値を表すパラメータεを特別に処理する機能を設計・実装した.徵小値の処理を含んだアルゴリズムは図３のようになる。通常のHydLaモデルのシミュレーションアルゴリズム[4]に加えて，高次項の削減( ReduceEp- silon (2.1節))，高次項の削減により起こる問題の解決(checkMinT, ShiftTme (2.2節))，不適切な場合分けの削減( GetElement (2.3節))，差分値の計算( DiffEpsilon (2.4節))，極限値の計算( LimitEpsilon (2.5節))が追加されている。以降でこれらの処理について説明する.
2.1高次項の削減
微小誤差を扱う際に，誤差が十分に小さいと考えると，εの2次以降の項に関しては興味がなく，無視しRequire:  HydLaブログラムHydLaProgram,シミュレーション終了時刻MaxT MS := TopologicalSort(SolveCH(HydLaProgram)) Mall := MaxModuleSet(MS) V := GetVariables(HydLaProgram) T := 0, S := ture, CP := ture while T <cp  MazT do // PP $ := SubstituteMinTime(S,T) ( S, CP, E, -, -)  := CalcMCS(S, MS, E, CP, T, CheckConsistencyPP) if S = false then break end if ( S, CP ) := AddParameters( S, CP, V ) S := ReduceEpsilon(S) // IP (S, CP, E, A, A_) := CalcMCS(S, MS, E, CP, T, CheckConsistencyIP) $ := SolveDifferentialEquation(S) if S = false V ¬IsUnique(S, V) then break end if S := ReduceEpsilon(S) getMinT := false , tmpMinT := 0 while !getMinT do ( MinT, CP ) := GetElement( CompMinTime ( { FindMinTime(S ∧ CP ⇒ g)|(g ⇒ c) ∈ A_ } ∪{ FindMinTime( S ∧ CP ⇒ ¬g)|(g ⇒ c) ∈ A + } ∪{ FindMinTime( S ∧ CP ∧ M_ ) | M_ ∈ ( Ma11 \ M ) } ∪{ FindMinTime( S ∧ CP ∧ ¬M÷)| M÷ ∈ M} U{(MaxT — T, true)})) getMinT := checkMinT(MinT,CP) (S, CP) := ShiftTime(S, CP, MinT) tmpMinT :=  tmpMinT + MinT end while MinT := tmpMinT T :=  MinT + T end while DiffS := DiffEpsilon(S) , S := LimitEpsilon(S)
図3　微小な値を扱うHyLaGIアルゴリズム
ても良い場合がある。このような場合のためにεに関するn火以上の高次項を削減する機能がある．この機能は，HyLaGIが数式処理を用いてシミュレーションを行っているため，複雑な式の計算などは処理が重くなり計算できなくなる場合もあるので，興味のない高次項を削減することで計算を軽くするという目的がある.まず，各フェーズにおいてそれぞれの変数はεを含んだ式f(ε)で表される。このf(ε)中に現れるεの高次項を削減するために，ここではε = 0の近傍でのn次近似を用いることとする。これは，f(ε)をε = 0においてテイラー展開したのちに，n + 1次以降の高次項を切り捨てて近似する方法である．ここでReduceEpsilonの処理はn > 0の時のみ適用する,
2.2高次項の削減による問題点の解決
高次項の削減により，各変数値は数式処理で求めた厳密な値からずれることになる．それによって，円状の枠の内部を運動していたボールが円の外側に出てしまったり，地面を跳ねていたボールが床を貫通してしまったりすることが考えられる。これらのような近似しない数式処理シミュレーションの結果と定性的に異なるふるまいを，望ましくないふるまいとして適切に処理する必要がある.例えば，円状の枠の内部を運動するポールの例題に関して考えるとプログラムは図４のようになる．ここで枠に衝突した瞬間を近似した際に，ボールの位置がずれることを考える。すると，枠の内部にずれたときは，近似しなかった場合と同様に枠で跳ね返ったのと同様の運動をするが，枠の外側はすぐにもう一度枠にぶつかって外側に跳ね返ってしまう.これは円状の枠の内部を運動するボールのシミュレーションとしては望ましくない.
図4  円の内部を跳ねる質点
本研究では次の離散変化時刻と成立するガード条件を用いて対処する手法を考案した．この手法は，(i)衝突後の処理において次の離散変化時刻を求めた際に，その時刻の極限ε → 0を求めると0になり，かつ(ii)その際の離散変化の原因となるガード条件が,衝突の原因となったガード条件と等しい場合，これを無視する方法である.近似により高次項が削減されて起こる各変数値の変化は，切り捨てられるεのn + 1次以降の項からなるため，極限ε → 0を考えると，ずれは限りなく小さくなっていくことになる。したがって(i)(ii)の場合に離散変化を無視する理由は，近似による変化によって，定性的に異なるふるまいを起こしてしまうからとなる.この方法では，三体衝突のような特殊点の近傍の解析の場合や1/εのような問接的に無限大を扱っている問題の場合も，望ましいふるまいかどうかを適切に判定できる．本研究ではこちらの手法を用いて対処した.本研究の過程で，別の手法として衝突後の状態とガード条件との関係を用いて動作を検査する手法についても考えた．しかし，この手法では望ましい運動との比較が必要なので，どうなると望ましい運動なのかあらかじめわかっておく必要がある。この円状の枠の内部を運動する例題だと，誤差ε = 0の場合が望ましい場合の運動だと考えられる。ところが，三体同時衡笑のような特異点の近傍における運動の解析をする場合，そもそもε = 0が特異点なので，単純にε = 0を代入した状況でガード条件との関係を比較すると,ボール同士がめり込んだりするなど，意図しない動作を起こしてしまう．また，1/εのような問接的に無限大を扱っている問題の場合，そもそもε = 0の場合を計算できないことがある。以上のように、ε＝０の場合が必ずしも望ましいふるまいになるとは限らず,望ましいふるまいを正確に求めることができない問適点があるため，本研究では採用しなかった．
2.3不要な場合分けの削減
徵小誤差εは，HyLaGI内では記号パラメータ“eps”として扱われている．HyLaGIは記号パラメータの値によって，動作の定性的な違いがある場合，場合分けを行ってシミュレーションする。微小誤差εも同様に場合分けが行われるが、εは十分に小さいものとして仮定しているため，本手法では，εの範囲が0もしくはその近傍を含まない場合の削減を行っている.
2.4巻分値の計算
徵小誤差が与える影響を考える時に、運動が進むにつれて誤差の影響がどのように拡大していくのか,または収束していくのかは，システムの性質を解析する上で重要な問題である．そこで，微小誤差を含んだ変数値の変化だけではなく，初期時刻における誤差がどう影響していくのかは，大変重要である．ここでは，εを含む元の式とε = 0を代入した式,つまり誤差εが0だった場合との差分値の計算を行う．この処理はシミュレーションの最後に行われるため，場合分けなどの問題は発生しない．
2.5極限値の計算
微小誤差が限りなく0に近づいた場合，つまり極限ε → 0において各変数がどうなるのかも，システムの重要な性質である。そこで，各フェーズにおいてそれぞれの変数値の極限ε → 0を求めた値も出力している．ここで1/εなどは無限大となるため，問接的に無限大を扱ったシミュレーションが可能となる.例えば，図５のようにすると撃力を微小値を用いて表現できる。これは静止しているx1にε時間だけ1/εの加速度が加わり，その後等速運動をするプログラムとなる.このようなプログラムを実行し,シミュレーションの後でこの極限ε → 0をとることにより,瞬間的に無限大の力が加わったシミュレーションがでめる.
3例題による評価
3. 1円の内部を跳ねる質点
ここでは，高次項の削減の説明(2.2節)で利用した図４のプログラムに関して，実行にかかる計算時間がどう変化するかを調べた。ここでこのプログラムの解析において，誤差εは十分に小さく，その2次以降は無視して良いという仮定のもとで，高次項を削減することを考える.まず高次項の削減を行うことにより，ボールが円周に衝突した瞬間の位置と時刻において，２次以降の項がもたらしていた微小な違いは削減されてしまう.そのため，各変数は微小にずれることになる．しかし，高次項(2次以降)を削減しても誤差の1次の項が残っているため誤差がどのような影響を与えているかは解析することができる.ここで計測したのは，実行時間の変化となる．図６は高次項の削減を行わない場合の各フェーズの計算にかかる時間で，図７は高次項の削減を行った場合の各フェーズの計算にかかる時間となる.この図6を見ると，高次項を削減しない場合ポールが円の枠に衝突するたびにだんだんと数式が複雑化するため，徐々に計算に時間がかかるようになっていくことがわかる。しかし，図７では高次項の削減を行っているため，εの2次以降の項が削減されるので，数式の複雑化が起こらないようになっている．そのため，各フェーズにかかる時刻は，ある一定以上には増えていないことがわかる.
3.2玉突き衝突
これはxy平面上のy軸上に距離1の問隔で半径1のポールが並んだモデルである。そこにポールをぶつけて玉突き衝突を起こしていく．ぶつけるボールはx軸方向に微小にずれていて，その誤差の影響がどのように伝達していくかを見ていく、プログラムは図8のようになる.このプログラムに対しHyLaGIの従来のシミュレーションを行うと，急激に数式が複雑化してしまい，4つ目のボールに衝突するまでのシミュレーションができないプログラムとなる。しかし，高次項を削減することで図9のようにシミュレーションを行うことがでめる.図9では，一つ目のボールが右に少しずれた位置から運動し始めるため２つめのポールは左，３つめのボールは右，のように左右交互にずれが発生するふるまいになることがわかる.
3.3ニュートンのゆりかご
これはx軸上を一直線に半径0.5のポールが4つ隙問なく並んでいて，そこに5つ目のボールをぶつけていく例題となる。ボールは等速直線運動を行い，衝笑は完全弾性衝突である．両端には壁があり，壁の間を5つのボールが運動するようなプログラムとなる.通常の方法ではn体同時衝突(n > 2)がおこり，シミュレーションができないプログラムになるが，図10のようにボールの大きさに微小値を加えることで, εの大きさの隙間を開けて並んでいるように近似すると，二体衝突が連続で起きるようになり，シミュレーションを行えるようになる.このプログラムのシミュレーションを行うと，図11のようになり，微小ながら時間差で順にポールが衝突していくことがわかる．従って，隙間なく並んだ
図10  ニュートンのゆりかごのHydLaプログラム
ボール４つにボールをぶつけていくと，反対側のボールに運動が伝わるという，ニュートンのゆりかごのふるまいをシミュレーションできた.また，微小値の極限ε ＝ 0を考えると図12のように理想的なニュートンのゆりかごも扱うことができた.
3.4 流される質点の運動
このモデルでは二次元空間での質点の運動を考える．図13のように質点が置かれている空間を，y軸と並行な3つの領域A, B, Cに分割し，それぞれの領域では質点に対し異なる力が働くものとする．領域Aではx軸方向に0.5,  y軸方向に−0.2,領域Bではx軸方向に0.01, y軸方向に−0.1,領域Cではx軸方向に−1.5, y軸方向に0.3の加速度で力が働いていると考える。質点の初期位置は点(5 + ε,  −10)とし，初速度は0とする．プログラムは図14のようになる. 2 3 4 5 6 7 8 9 10この質点の運動の軌道は次の図15のようになる.この軌道で初期値の誤差がどのような影響を与えるか分析する．微小値の与える誤差の差分値を計算すると図16のようになる。図16では，各時刻におけるxとyの本来の軌道(誤差ε=0の場合)とのずれを表している．赤線がｘで青線がｙを表している．これはεが1まで変化した場合の軌道のずれを描いてある．誤差の影響が小さくなっている部分を問ってあるが，特定の時刻において影響が小さくなっていることがわかる.このような流される質点の運動において，質点の運動に関してある特定の時刻で初期値の誤差εの影響が小さくなることが分析できた.このように誤差の影響が収束するような点が分析できるため，物体の運動に不確定性がある場合に，システムの性質を利用した応用が考えられる。例えば,ロポットに物体をキャッチさせたい場合に，物体の運動の不確定性が収束するようなキャッチ可能な地点を計算する，などの応用が可能になる.
4  関連研究
文献[2]では，従来のHydLa，HyLaGIで直接扱うことのできない離散変化が連続しておこるモデルを扱っている。この文献中の例題Fig.8,9 (ニュートンのゆりかご)やFig.11,12 (ポールが左右からほぼ同時にぶつかるニュートンゆりかご), Fig.14 (ボールの質量がすべて異なるニュートンのゆりかご)は，本論文の図10と同様に微小値を用いたモデリングを行い，解析することができた．また，Fig.15,16 (２つ並んだポールを左から押す例題)，Fig.28 (上下に垂なった物体の下の物体に撃力を与える例題)は，微小値を用いることでモデルをHydLaプログラムで表すことができたが，HyLaGIの実装が完全ではないため，まだ解析は行えていない．原理的にはシミュレーション可能な例題であるため，HyLaGIを修正でき次第解析を行う予定である.
5　まとめと考察
5.1考察
提案手法を用いることで，微小値であることを利用して高次項の削減を行えるため，だんだんと数式が複雑化していくシミュレーションにおいてその複雑化を軽減して計算時間の削減を行うことができた。また,同様に数式が複雑化することによって，通常だとシミュレーション不可能な例題もシミュレーションを行うことができた．高次項の削減は，興味のない微小値のn次以上の項を削減するだけでなく，数式を簡単にすることで，HyLaGIのパフォーマンスの向上という効果もあった.また，ニュートンのゆりかごのような特異点の解析の際に利用可能で，微小値が限りなく0に近いということを利用して，特異点近傍での運動を考察できた．同様に無限大や無限小を扱うような表現しにくいシステムもモデリング可能であると考えられる.さらに，流される質点の運動のシミュレーションのように，初期値の誤差が軌道に与える影響を調べることができ，誤差の収束や拡大などシステムのもつ性質を分析した応用に利用できる.
5.2まとめ
本論文では，HyLaGIにおいて微小値としてパラメータεを特別に処理する機能の評価を行った．微小な値を扱う機能がHyLaGIのシミュレーション能力を向上できたことや実際にモデルに誤差が与える影響から，システムの解析に利用できることがわかった.今後の課題として，４章の文献[2]中の例題で，まだ解析できていない例題があるため，これを解析していきたいと思う。また，複数の微小値を用いた解析のために機能拡張を行い，その際にシステムの性質の解析や考察を容易にするために，それぞれの微小値について適切に出力する方法など考察していきたい.謝辞本研究を行うにあたり，例題作成に協力していただいた皆様，貴重な意見をいただいた上田研究室HydLa班の皆様に感謝する。本研究の一部は，科学研究費補助金(基盤研究(B) 26280024)の補助を得て行った.