2018 年度 修士論文
非線形な常微分方程式を含む ハイブリッドシステムの パラメータ付き精度保証計算
提出日 : 2019 年 2 月 1 日 指導 : 上田 和紀 教授
早稲田大学 基幹理工学研究科 情報理工・情報通信専攻
学籍番号 : 5117F081-7
増田 健太
概要
ハイブリッドシステムとは，複数の状態変数と常微分方程式によるそれらの連続変化，
そして離散変化により構成される動的システムのことである．これはセンサやアクチュ エータを備えて物理環境と相互作用を起こすサイバーフィジカルシステムのようなシステ ムのモデル化に用いることができる．HyLaGI はハイブリッドシステムモデリング言語
HydLa の C++ による処理系で，Mathematica を用いた記号的解析を行うという特徴
を持っている．HyLaGIはその特徴から，ハイブリッドシステムの厳密な解軌道を求める ことができ，また初期値にパラメータを含むようなモデルについてもそのパラメータを記 号として残せるという利点を持つ一方で，解析的に解けない非線形な常微分方程式を扱う ことができないという欠点も持つ．本研究は，非線形な常微分方程式を含むハイブリッド システムを，HyLaGIで解ける形に過大近似されたハイブリッドシステムに変換し，精度 保証解を求められるようにすることを目的とする．過大近似を行うため HyLaGI のハイ ブリッドシステムの厳密な解軌道を求めるという特徴は失われるが，初期値に含まれるパ ラメータを記号として残すという特徴は保たれるので，解の初期値に対する依存性を調べ ることはできる．
提案手法ではまず常微分方程式を構成する非線形関数について，平行体の集合を用い てこれを包含することで，幅を持つ線形微分方程式の集合に変換する．次に，この集合 を平行体同士の遷移を離散変化として表現したハイブリッドシステムに変換し，これを
HyLaGI の区間モードを用いることで解析を行う．以上の手法を非線形モードとして
HyLaGIに実装し，例題を用いて実験及び評価を行った．実験の結果として，非線形関数
が1変数および2変数からなる常微分方程式について，初期値に含むパラメータを記号と して保存した精度保証解を計算することに成功した．
Abstract
A hybrid system is a dynamical system composed of state variables, their contin- uous dynamics by ODEs (Ordinary Differential Equations) and discrete transitions.
This can be used for modeling systems such as cyber-physical systems equipped with sensors and actuators and interacting with the physical environment. HyLaGI is an interpreter of the hybrid system modeling language HydLa and features symbolic anal- ysis using Mathematica. HyLaGI can simulate exact behaviour of hybrid systems and can preserve parameters in initial values as symbols in solutions. However, HyLaGI cannot deal with nonlinear ODEs which cannot be analytically solved. The purpose of this research is to calculate the guaranteed-accuracy solution to a hybrid system containing nonlinear ODEs by converting it to an over-approximated hybrid system.
Although we cannot obtain an exact solution to a hybrid system because of over- approximation, parameters in the initial values of the hybrid system are preserved as symbols in the solution. Therefore, it is possible to investigate the dependence between the solution and the initial value.
In the proposed method, first we obtain the set of linear ODEs by calculating parallelotopes enclosing the nonlinear function constituting the original ODEs. Next, we convert these linear ODEs to a hybrid system that has discrete changes that represent transitions between each parallelotope. Then, we analyze this hybrid system using HyLaGI’s interval mode. We implemented the above method in HyLaGI as a nonlinear mode, and we experimented and evaluated it using examples. As a result of the experiment, for the ODEs consisting of one variable and two variables for nonlinear functions, we succeeded in computing the guaranteed-accuracy solutions which preserved parameters in initial values as symbols.
i
目次
第1章 はじめに 1
1.1 研究の背景と目的 . . . 1
1.2 論文の構成 . . . 2
第2章 ハイブリッドシステムとその解析手法 3 2.1 ハイブリッドシステム . . . 3
2.2 記号シミュレータ HyLaGI . . . 4
2.3 精度保証シミュレータ Flow* . . . 5
第3章 非線形な常微分方程式の平行体を用いた変換と解析 6 3.1 提案手法の概要 . . . 6
3.2 モデル中の非線形常微分方程式の検出 . . . 6
3.3 平行体による非線形関数の包含 . . . 8
3.4 ハイブリッドシステムへの変換と置き換え . . . 10
3.5 HyLaGI での実行 . . . 11
3.6 区間モードの改善 . . . 12
第4章 実験と評価 14 4.1 初期値が固定値のみからなる例題 . . . 14
4.2 初期値にパラメータを含む例題 . . . 17
第5章 関連研究 20 第6章 まとめと今後の課題 21 6.1 まとめ . . . 21
6.2 今後の課題 . . . 21
目次 ii
謝辞 23
参考文献 24
発表論文 25
iii
図目次
3.1 提案手法のフローチャート . . . 7
3.2 2変数の非線形関数 . . . 9
3.3 四分木による空間分割 . . . 9
3.4 平均値形式による評価の誤差 . . . 9
3.5 十分狭い幅を持つ式の求根 . . . 13
3.6 広い幅を持つ式の求根 . . . 13
4.1 単振り子（パラメータなし）の実験結果 . . . 15
4.2 ファンデルポールの方程式（パラメータなし）の実験結果 . . . 16
4.3 単振り子（パラメータあり）の 1/4周期前後における解軌道 . . . 18
4.4 単振り子（パラメータあり）の半周期前後における解軌道 . . . 18
1
第 1 章
はじめに
1.1 研究の背景と目的
ハイブリッドシステムとは，複数の状態変数と常微分方程式によるそれらの連続変化，
そして離散変化により構成される動的システムのことである．HyLaGI はハイブリッド システムモデリング言語 HydLa のC++ による処理系で，Mathematica を用いた記号 的解析を行うという特徴を持っている．HyLaGIはその特徴から，ハイブリッドシステム の厳密な解軌道を求めることができ，また初期値にパラメータを含むようなモデルについ てもそのパラメータを記号として残せるという利点を持つ一方で，解析的に解けない非線 形な常微分方程式を扱うことができないという欠点も持っている．
本研究は非線形な常微分方程式を含むハイブリッドシステムを，HyLaGIで解ける形に 過大近似されたハイブリッドシステムに変換し，精度保証解を求められるようにすること を目的とする．過大近似を行うので HyLaGI のハイブリッドシステムの厳密な解軌道を 求めるという特徴は失われるが，初期値に含まれるパラメータを記号として残すという特 徴は保つことができ，解の初期値に対する依存性を記号的に調べることができる．
非線形な常微分方程式を含むハイブリッドシステムの精度保証された解軌道を求める先 行研究としては Flow* が挙げられる．Flow* では常微分方程式をテイラーモデル法とい う多変数の多次元多項式を用いた解析手法を用いることにより，解の高精度な包含軌道を 計算している．しかし，テイラーモデル法では浮動小数点数の計算誤差と，初期値に含ま れるパラメータの誤差をまとめて一つの誤差幅として扱うことから，パラメータを記号と して残すことはできないという問題がある．特に，複数のパラメータを持つシステムの場 合，解がそれぞれのパラメータにどの程度影響されるかを調べるのはこのような数値計算 に基づく精度保証計算では困難である．
第1章 はじめに 2
また，HyLaGI における先行研究として，離散変化時刻が記号的に求められないよう
なハイブリッドシステムについて，区間ニュートン法とアフィン演算を用いることでパラ メータ間の一次の依存関係を保存した精度保証解を計算する区間モードという手法が存在 する．この手法では精度保証計算を用いているものの，常微分方程式は解析的に解いてい るので，非線形な常微分方程式を扱うことはできない．
提案手法ではまず非線形な常微分方程式について，それを構成する非線形関数を複数の 平行体で包含することにより，定数項に区間を持つような線形関数の組に変換する．ここ で，常微分方程式の非線形関数をその線形関数で置き換えた線形微分方程式は，それぞれ の定義域において元となる非線形常微分方程式の過大近似となっている．そして，これら の線形常微分方程式を連続変化として持ち，それらの切り替えを離散変化で表現したハイ ブリッドシステムを構成する．これを記号的に解析することにより，パラメータを残して 精度保証された解軌道を計算することが可能となる．パラメータが記号として残すことが できれば，個々のパラメータが解に与える影響と，計算誤差による影響をそれぞれ分離し て扱うことができ，解の初期値に対する依存性を解析することが容易になる．また，パラ メータの制御問題への応用も考えられる．
1.2 論文の構成
2章ではハイブリッドシステムとその解析ツールである HyLaGIと Flow* の説明を行 う．3章では非線形な常微分方程式を含むハイブリッドシステムを HyLaGI で解くため の提案手法について述べる．4章では提案手法の結果について，例題を用いて実験し，結 果と考察と評価を行う．5章では本研究と関連研究の関係について述べる．6章では本研 究のまとめ及び今後の課題について述べる．
3
第 2 章
ハイブリッドシステムとその解析 手法
本章ではハイブリッドシステムの説明とその解析ツールである，HyLaGIおよびFlow*
[6]についての説明を行う．
2.1 ハイブリッドシステム
ハイブリッドシステムとは，常微分方程式による連続変化と，ある条件が満たされた時 に発生する離散変化を含むような動的システムのことである．
ここでは，ハイブリッドシステムの解析手法を説明するにあたり，まずその解析の難し さにより大きく三種類に分類する．
2.1.1 線形ハイブリッドシステム
線形ハイブリッドシステムとは，全ての微分方程式の解が線形な式で表現できるような ハイブリッドシステムのことである．このようなハイブリッドシステムについては，記号 的に到達可能性解析を行うことができる．HyLaGIでは，この種類のハイブリッドシステ ムについて，記号解析を行うことが可能である．
2.1.2 非線形ハイブリッドシステム
非線形ハイブリッドシステムとは，微分方程式の解に非線形な式を含むようなハイブ リッドシステムのことである．このようなシステムについては，連続変化は記号的に解析
第2章 ハイブリッドシステムとその解析手法 4 できるものの離散変化の解析が記号的には行えない場合がある．HyLaGIでは，この種類 のハイブリッドシステムについて，区間モード[2]を用いることで精度保証付き解析を行 うことができる．
2.1.3 非線形な常微分方程式を含むハイブリッドシステム
常微分方程式自体が非線形なシステムは，連続変化も離散変化も一般には記号的に行う ことができない．Flow* では，テイラーモデル法[3]によりこのようなシステムの精度保 証付き解析を行うことができる．
2.2 記号シミュレータ HyLaGI
HyLaGI はハイブリッド制約言語 HydLa [1]の処理系であり，ハイブリッドシステム
の記号的解析を行うソフトウェアである．HyLaGI では連続変化となる常微分方程式が 線形でなおかつ離散変化時刻が記号的に解ける場合，そのハイブリッドシステムの厳密な 解軌道を求めることができる．HyLaGIの特徴として，初期値に区間を含むパラメータを 持つことができ，そしてそのパラメータを記号として残した解を求めることができるとい う点が挙げられる．
2.2.1 ハイブリッド制約言語 HydLa
HydLa とはハイブリッドシステムのモデリングのための制約宣言型言語である．
HydLa では連続変化と離散変化を全て制約として宣言することで記述し，それぞれの制
約の優先度を制約階層によって表現することで，制約を満たす解軌道の集合を定義する．
HydLa を用いたハイブリッドシステムのモデリング例を次に示す．
ソースコード2.1 HydLa モデルの例
1 INIT <=> y = 10 & y’ = 0.
2 FALL <=> [](y’’ = -10).
3 BOUNCE <=> [](y- = 0 => y’ = -4/5 * y’-).
4
5 INIT, FALL << BOUNCE.
第2章 ハイブリッドシステムとその解析手法 5
2.2.2 区間モード
常微分方程式は線形であるが離散変化時刻が記号的に解けないシステムについて，区間 モードを使用することによりハイブリッドシステムの精度保証された解軌道を数値的に求 めることができる．区間モードでは精度保証計算に kv ライブラリ[7]を使用する．
kv ライブラリ
kv ライブラリとは，C++ で実装された精度保証付き数値計算ライブラリのことであ る．kv ライブラリでは，べき級数演算[8]やピカールの逐次近似法を用いて，常微分方程 式の初期値問題の精度保証計算を行うことができる．HyLaGIでは，アフィン演算を用い て式の過大近似を行うのに kv ライブラリを使用している．
2.3 精度保証シミュレータ Flow*
Flow* はテイラーモデル法を用いてハイブリッドシステムの精度保証解析を行うソフ
トウェアである．
2.3.1 テイラーモデル法
テイラーモデル法とは，テイラーモデルを用いて非線形な常微分方程式の精度保証計算 を行う数値計算手法である．この手法はピカールの逐次近似法をテイラーモデルを用いて 拡張した数値計算手法で，精度保証計算特有の精度の低下を防ぐために，シュリンクラッ ピングやプリコンディションといった手法を用いた工夫を行っている[4][5]．テイラーモ デルとは，k 変数n次の多項式と区間を含んだ剰余項の和として表される式のことであ る．また，テイラーモデル法はkvライブラリにおけるピカールの逐次近似とは異なり，
ステップをまたぐときに解を単一の区間値に評価するのではなく，記号を残した状態で次 のステップの初期値にするという特徴がある．
6
第 3 章
非線形な常微分方程式の平行体を用 いた変換と解析
本章では，本論文の提案手法について説明する．
3.1 提案手法の概要
提案手法は，HydLa モデルとして記述されたハイブリッドシステムを入力として，モ デル中に含まれる非線形常微分方程式を，複数の幅を持った線形常微分方程式で包含して ハイブリッドシステムに変換し，それを HyLaGI によって解析するものである（これを 以降では非線形モードと呼ぶ）．なお，本手法では非線形常微分方程式は1変数または2 変数のみからなるものとする．提案手法の概要は図 3.1 に示す．手順は主に，HydLa モ デル中に含まれる非線形常微分方程式の検出と，平行体を用いた非線形関数の包含，包含 した平行体のハイブリッドシステムへの置き換え，そして HyLaGI での実行の4ステッ プからなる．以下にそれぞれのステップと既存手法に加えた改善点の説明を行う．
3.2 モデル中の非線形常微分方程式の検出
まず，非線形常微分方程式を含む HydLa モデルの記述の仕方についての説明を行う．
HydLa 言語では微分方程式の記述の仕方として，変数の微分を含む任意の方程式を構文
として認めており，実際に HyLaGI でも多くの微分代数方程式を解くことができる．し かし，非線形モードでは与えられた微分方程式が非線形であることを検出するために，微 分方程式の形について，等式の左辺が変数の微分のみで構成され，右辺が微分を含まない
第3章 非線形な常微分方程式の平行体を用いた変換と解析 7
図3.1 提案手法のフローチャート
式のみで構成される x′ =f(x)のような形式に制限する（なお，線形常微分方程式と非線 形常微分方程式の両方を含むモデルを扱うとき，線形常微分方程式の記述はこの形式に従 う必要はない）．
ソースコード3.1 非線形モードの例
1 INIT <=> x = 0.
2 UPDATE <=> [](x’ = x^2 + 1).
3
4 INIT, UPDATE.
5 // #hylagi --fnonlinear_vars x
次にこの入力の HyLaGI における解釈の仕方について説明をする．HyLaGI でこの コードを受け取るとると，まず –fnonlinear vars オプションにより非線形モードが有効 になる．そして，パース直後に構文木を走査し，左辺が微分値で右辺が –fnonlinear vars に指定された変数についての非線形な式であった場合，この微分方程式を非線形モードに 登録し，構文木を線形な式に置き換える．
第3章 非線形な常微分方程式の平行体を用いた変換と解析 8
3.3 平行体による非線形関数の包含
次に，常微分方程式を構成する非線形関数を，複数の平行体によって包含する方法につ いて説明をする．平行体とは，対辺が平行な N 次元の立体のことで，二次元であれば平 行四辺形，三次元であれば平行六面体となる．これは，ある定義域の変数値に対して，非 線形関数を平均値形式で評価することにより得られる．ただし，誤差を少なくするために は平行体の幅をできるだけ狭くする必要があるため，空間を木構造に分割し分割位置や包 含幅の最適化を行う．この手順について順を追って説明をする．
3.3.1 平均値形式による平行体の計算
まず常微分方程式を構成する非線形関数について，この関数を平均値形式の形でアフィ ン演算を用いて評価する．
平均値形式とは，次のようにして定義される．
fm(x) :=f(xm) +f′(x)(x−xm) (xm ∈x) (3.1) ここで，x は区間ベクトル， xm はx に含まれるベクトル，f は区間ベクトルから区 間ベクトルへの写像である．ただし，包含を行う平行体が一つだけだと幅が広すぎて誤差 が実用的な範囲に収まらないため，空間を分割して小さい区間ごとに平行体で包含するこ とを次に考える．
3.3.2 木構造による空間分割
空間分割は，分割位置の最適化と隣接要素の列挙を容易にするために，非線形関数が1 変数関数ならば二分木，2変数ならば四分木を用いて再帰的に分割を行う．
例として，図 3.2 で示す2変数の非線形関数を分割すると図 3.3 のようになる．
3.3.3 微細な包含の計算
上記の空間分割を繰り返すことで，任意の精度の平行体の集合を作ることができる．し かしながら，平行体同士の間を遷移する離散変化時刻は区間ニュートン法によって求める ため，分割数を増やすほど浮動小数点計算による誤差が累積してしまう．したがって，精 度の良い結果を得るためにはできるだけ少ない分割数で幅の狭い平行体の集合を作る必要 がある．
第3章 非線形な常微分方程式の平行体を用いた変換と解析 9
図3.2 2変数の非線形関数 図3.3 四分木による空間分割
図3.4 平均値形式による評価の誤差
ここで，平均値形式の評価により発生する誤差について考える．図 3.4 に示す外側の平 行四辺形は，1変数関数を二つの区間に区切って，それぞれの区間で評価したものである．
これを見ると，上下に隙間があり明らかに最適な幅ではないことがわかる．これは平均値 形式での評価というのは，直感的には区間中の最も傾きが小さい方向と最も傾きが大きい 方向に広げて交差する領域を取るような計算なので，区間中の傾きの変化が大きいほど余 分な幅を持ってしまう．
そこで，平行体を計算した後に余分な幅を削ることを考える．これは，平行体と非線形 関数の交差判定が行えれば可能だが，1変数である場合を除いて非線形関数のままで交差 判定を行うのは難しい．したがって，計算過程を二回に分け，一度目では微小幅まで空間 分割を繰り返し，幅がゼロに近い平行体の集合を得る．そして，二度目の計算で平行体を 計算した後，非線形関数の代わりに一度目に計算した平行体の集合を用いて交差判定を行
第3章 非線形な常微分方程式の平行体を用いた変換と解析 10 うことで，幅の削減を行う．
3.3.4 平行体の包含幅の最適化
前のステップで得られた微細な平行体の集合を用いて，平行体の幅を最適化することを 考える．これは，傾きを決めてしまえば，切片は平行体の集合に接する値が決まるので，
最適な傾きを求められれば良い．本手法では，すべての平行体の頂点についての最小二乗 平面を求め，この傾きを採用することにした．
3.3.5 分割位置の最適化
本手法では関数が2変数ならば四分木を用いて空間分割を行うが，この時中央で分割す るのが最適とは限らない．したがって，最適化計算により最も平行体の体積の合計が小さ くなるような分割点を求め，この点で分割を行う．分割点の選び方に対して体積の変化の 仕方は穏やかであるため，最適化手法には収束の速い準ニュートン法を用いた．
3.4 ハイブリッドシステムへの変換と置き換え
ここまでで，非線形関数を包含する線形関数の集合が得られた．あとは，それぞれの平 行体同士についてどこからどこへの遷移が発生し得るか，つまり，どの平行体同士が隣接 しているかを求められればハイブリッドシステムに変換することができる．
3.4.1 隣接要素の列挙
隣接要素の列挙は，分割位置の最適化を行っていなければ，空間分割を行った木の構造 情報のみを用いて行うことができる．しかし，それぞれのノードで異なる分割位置を取る 場合は，構造的には接していないが実際には接しているような場合が考えられる．した がって，構造情報から接している可能性のある辺同士をまず列挙し，辺同士について区間 の比較を行って実際に重なる領域があるかを判定することで隣接判定を行う．
3.4.2 ハイブリッドシステムへの変換
以上の情報を用いて，非線形な常微分方程式を幅を持った線形な常微分方程式からなる ハイブリッドシステムに変換することができる．これには，平行体を識別するための離散 的な状態変数を一つ導入し，その値についてのガード付き制約として微分方程式を記述
第3章 非線形な常微分方程式の平行体を用いた変換と解析 11 し，平行体間の遷移を状態変数の離散変化として記述する．
これによって変換したハイブリッドシステムの例を次に示す．
ソースコード3.2 線形化したモデルの例
1 INIT <=> x = 0.
2
3 Update0 <=> [](p- = 0 => x’ = x*A0 + B0 + C0*w).
4 Update1 <=> [](p- = 1 => x’ = x*A1 + B1 + C1*w).
5 ...
6
7 JumpL1 <=> [](p- = 1 & x- = S0 => p = 0).
8 JumpL2 <=> [](p- = 2 & x- = S1 => p = 1).
9 ...
10
11 INIT, 0 <= w <= 1, [](w’ = 0), ([]p’ = 0, Update0, Update1, ...) << (JumpL1, JumpL2, ...).
3.5 HyLaGI での実行
上記の手法の結果を HydLa モデルとして出力すると，HyLaGIの非線形モードではな い通常の区間モードで実行することが可能である．しかし，この方法ではごく簡単なモデ ルについてしかうまく解を求めることができないという問題がある．この原因は主に，パ ラメータによる区間ニュートン法のエラー，パラメータによる分岐発生時のエラー，変換
した HydLa モデルの実行速度の問題，の三つが挙げられる．したがって，非線形モード
では上記の結果を HydLa モデルには変換せずに HyLaGI で直接扱うアルゴリズムを組 み込むことで対処する．
3.5.1 パラメータによる区間ニュートン法のエラー
離散変化時刻を区間ニュートン法で求めるとき，結果の正しさを保証する条件として解 の唯一性の証明を行う必要がある．これは解軌道の区間幅が広いほど失敗しやすいため，
現実的な大きさのパラメータではそのまま実行することができない．したがって，失敗が 起きる場合はパラメータを分割することで対処を行う．ただし，パラメータの分割は誤差
第3章 非線形な常微分方程式の平行体を用いた変換と解析 12 が累積してからでは遅いので，シミュレーションを始める前に必要な数だけ分割を行うこ とにする．このため，非線形モードの実行では解の唯一性証明に失敗するたびに，パラ メータの大きさを半分に分割してシミュレーションを最初から実行し直す．
3.5.2 変換した HydLa モデルの実行速度の問題
本手法の結果をそのまま HydLa モデルとして出力すると，あるガードが満たされる まで考慮する必要のない制約や，解く必要のないただの代入式などについて全て解こう としてしまうため，非常に実行速度が遅くなってしまうという問題がある．したがって，
HyLaGI内部でハイブリッドシステムを作成した後，微分方程式中の平行体の値を実行中
に動的に書き換える．また，離散変化についても同様に，現在の状態変数の値を見て満た し得るガード条件を持つモジュール制約を動的に追加する．
3.6 区間モードの改善
区間モードでは区間ニュートン法を用いて離散変化時刻を導出する．区間ニュートン法 の精度を調べるために幅を持ったパラメータa, bについてf(x) =−(sin(x)−ax+b)の 根を求める実験を行った．
まず，a ∈ [−0.23,−0.05], b∈ [0.77,0.85]とした時の区間ニュートン法の結果を図 3.5 に示す．これを見ると結果として十分に狭い範囲の区間が得られていることがわかる．
次に，a ∈ [−0.23,−0.05], b∈ [0.69,0.93]とした時の区間ニュートン法の結果を図 3.6 に示す．こちらは先ほどとは異なり，区間ニュートン法の解として得られる範囲が必要以 上に広くなってしまっている．
以上の予備実験により，広いパラメータを持つ式について区間ニュートン法は解として 過剰に広い結果を返す場合があることがわかった．
区間ニュートン法は解の唯一性保証を行うため，この結果はただ一つの連続した区間で あることがわかっている．したがって，両端からゼロ点を含まない範囲を二分探索によっ て削ることで，精度を保証したままこの結果を改善することが可能である．この改善を実 装した結果，図 3.6 の黄緑の区間を得ることができた．
第3章 非線形な常微分方程式の平行体を用いた変換と解析 13
図3.5 十分狭い幅を持つ式の求根
図3.6 広い幅を持つ式の求根
14
第 4 章
実験と評価
本章では，提案手法について例題を用いて実験を行い評価する．実験は，非線形な常 微分方程式を含む問題について，初期値にパラメータを含む場合と含まない場合でそれ ぞれ実験を行った．実験に使用した環境は，CPU：Intel(R) Xeon(R) CPU E7-4860 @ 2.27GHz，メモリ：256GB である．
4.1 初期値が固定値のみからなる例題
今回実験に用いたハイブリッドシステムは，動的にパラメータが導入されるような変化 を含まないので，初期値が全て固定値であれば，その厳密解は時刻に対して一意に定ま る．ただし，解析には精度保証計算を用いるため，実際の解軌道は時刻に対して区間を持 つ値となる．したがって，この区間幅が誤差となる．
4.1.1 単振り子の運動
1変数の非線形微分方程式の例題として，単振り子の運動方程式について実験を行った．
d2θ(t)
dt2 =−g
l sinθ(t) (4.1)
dθ(t) dt
t=0
= 0 (4.2)
実験は初期角度θ(0) = 20[deg]の場合とθ(0) = 60[deg]の場合についてそれぞれ空間 分割の数を変えながら，振り子の周期を求めた．定数値はそれぞれg = 9.8，l = 5 と した．
結果を次に示す．
第4章 実験と評価 15
図4.1 単振り子（パラメータなし）の実験結果
グラフの横軸が分割に使用した平行四辺形の数，縦軸が周期の区間幅の狭さを表してい る．グラフより，θ(0) = 20[deg]とθ(0) = 60[deg]では常に 2桁程度精度が異なり，分 割数を増やした時の精度の上がり方はおおよそ同じ傾向を示すことが分かった．また，近 似曲線を見ると，分割数を2倍〜3倍にすると精度が一桁程度向上するということがわ かった．
4.1.2 ファンデルポールの方程式
2 変数の非線形微分方程式の例題として，ファンデルポールの方程式について実験を 行った．
dq
dt =p (4.3)
dp
dt =K(1−q2)p−q (4.4)
初期値はp(0) = 0, q(0) = 1，定数値はK = 0.25として実験を行った．
結果として得られた解軌道を図 4.2に示す．
分割はp, qともに[−3,3]の範囲で分割数300000に設定して行ったところ，PP51（時
刻t= 2.18）までのシミュレーションに成功した．シミュレーションの失敗時には三つの
ケースに分岐していたことから，失敗の原因は区間幅の増大により解軌道が複数のガード
第4章 実験と評価 16
図4.2 ファンデルポールの方程式（パラメータなし）の実験結果
の境界をまたいでしまったことによるものと考えられる．
ファンデルポールの方程式には解が時間とともにリミットサイクルに吸い込まれるよう な振る舞いを起こし，次第に一定の周期運動に収束するという性質がある．したがって連 続軌道上では区間幅はむしろ縮小方向に向かうはずで，区間幅が拡大するとしたら離散変 化が原因となるはずである．
結果の軌道を見ると，連続軌道において区間が収縮しているかは確認できないが，区間 幅の拡大の仕方が段階的であることは確認できるので，実際に離散変化時刻の計算が原因 で精度が落ちているものと考えられる．
第4章 実験と評価 17
4.2 初期値にパラメータを含む例題
パラメータを含まない場合は区間幅が全て誤差だったので，単に解軌道の区間幅だけを 見て狭いほど良いという評価ができた．一方で初期値にパラメータを含む問題については 得られる結果の形式は変わらないものの，評価の仕方は異なる．結果の解軌道は誤差によ る幅とパラメータによる幅を持っており，誤差による幅は狭いほど良いがパラメータによ る幅は広くても問題ない．
4.2.1 単振り子の運動
前節で行った単振り子の運動方程式について，振り子の長さをパラメータとし，5 ≤ l ≤ 5.0625の範囲で振り子の軌道を求めた．初期角度は θ(0) = 10[deg]，重力加速度は g= 9.8，分割数は50とした．
実験の結果，振り子の半周期まで解を求めることができた．半周期時点での振り子の軌 道の式を以下に示す．
θ(t)∈[−0.16121,−0.16041] cos(angle(t))+
[−0.0060,−0.0054] sin(angle(t))+
[−0.966,−0.864]
√ 1 161 +len+ [0.001562,0.001568]
angle(t) =
√309.232225
161 +len ([−1.966,−1.935] +t−0.0061len)
(4.5)
ここで，lenという記号が振り子の長さを表すパラメータで，−1 ≤ len≤ 1の範囲を 取り5≤l ≤5.0625に対応する．この結果を見ると，解が振り子の長さlと時刻tについ ての関数として求められたことがわかる．この式から読み取れることは，θ(t)の後ろの2 つの項は時刻tに依存していないので定数項で，その中でも3つ目が振り子の長さに依存 して変化する成分で，4つ目が純粋な誤差である．また，前半2つが時刻によって変化す る成分で，特にangle(t)の中身を見ると時刻tの係数である，
√309.232225
161+len の部分が角振
動数に相当する成分であることがわかる．
以上が軌道の式から読み取れることだが，実際に解がパラメータについての依存性をど れだけ記号として保存できたかを調べるために，解に具体的な数値を代入した結果をプ ロットした．図4.3 に1/4周期時点（振り子が中央を通る瞬間）での解軌道を，図 4.4 に 半周期時点（振り子が端にたどり着いた瞬間）での解軌道をそれぞれ示す．
第4章 実験と評価 18
図4.3 単振り子（パラメータあり）の1/4周期前後における解軌道
図4.4 単振り子（パラメータあり）の半周期前後における解軌道
それぞれの図を見ると 1/4 周期時点では，パラメータの下限と上限で全く異なる場所 を通っているため，誤差は十分に小さくパラメータの依存性を保存できていると言える．
しかし，半周期時点では区間幅が拡大しパラメータの下限と上限でほとんど区別がつかな いような結果になってしまった．さらに，この半周期前後の連続変化における区間幅の増 大の仕方は大きく，このフェーズの始点と終点では倍以上に幅が拡大していた．
この原因は主に二つ考えられる．一つは，振り子の加速度である sin(θ)は中央から遠 ざかるほど非線形な運動を起こすので，両端が線形に包含した時の幅が最も大きくなりや すいという点である．そしてもう一つは，両端で振り子の速度が0になるため，両端には
第4章 実験と評価 19 中央付近よりも長い時間滞在することになるという点である．この二点が原因となり両端 付近ではそもそも加速度の誤差が大きく，そしてそれが長い時間累積したため結果的に区 間幅が大きく増大してしまったのだと考えられる．
20
第 5 章
関連研究
本研究の関連研究として，ハイブリッドシステム解析器 Flow*および HyLaGIの区間 モードが挙げられる．
Flow* は非線形な常微分方程式を含むハイブリッドシステムの精度保証付き解析を行
うという点で本研究と共通している．しかし，Flow*は微分方程式の解の計算にテイラー モデル法を用いており，解には記号が残されていない．一方で，本手法では非線形な微分 方程式を線形な微分方程式に変換して記号的に解いているため，パラメータを記号として 残すことができている．
また，HyLaGI の区間モードは本研究の先行研究であり，非線形で記号的に解くこと
ができない問題について精度保証計算を用いて過大近似を行うという点で本研究と共通し ている．区間モードでは，記号的に解けない離散変化時刻について，区間ニュートン法を 用いることで真の最小離散変化時刻を含む時刻の区間を計算している．区間ニュートン法 の問題点は数値的な計算手法であるため，この結果として得られる区間が数値計算の誤差 とパラメータの影響のどちらも含んでしまうことにあるが，区間モードでは平均値の定理 に基づいて解からパラメータに関する線形な依存性を取り除くという処理を行っているた め，誤差とパラメータの影響の分離に成功している．とはいえ，パラメータの幅が大きく なると，区間ニュートン法で得られる区間の幅も広くなり，結果として区間ニュートン法 の結果が正しく精度保証されているか確認するための解の唯一性の保証に失敗してしま う．したがって，本手法では大きなパラメータについては事前に必要な数だけ分割してお き，それぞれの小区間について区間ニュートン法を用いて解を求めるという方法で計算を 行っている．
21
第 6 章
まとめと今後の課題
6.1 まとめ
本研究では，非線形な常微分方程式を含むハイブリッドシステムについて，入力に含ま れる非線形な常微分方程式を，幅を持った線形な常微分方程式から構成されるハイブリッ ドシステムに変換し，HyLaGI 上でシミュレーションする手法を提案及び実装した．ま た，実際に例題を用いて実験，評価を行った．結果として，非線形関数が1変数または2 変数からなる常微分方程式を含むハイブリッドシステムについて，パラメータに関する一 次の依存性を保存した精度保証解を計算することに成功した．しかしながら，誤差による 解軌道の変動幅に比べてパラメータによる解軌道の変動幅は小さく，計算誤差とパラメー タの記号的成分を十分に分離できたとは言えない結果であった．
6.2 今後の課題
今後の課題は主に三つ挙げられる．それぞれについて説明を行う．
6.2.1 ガードの境界をまたぐ場合の考慮
広いパラメータを扱うと解軌道がガードの境界を包含するケースがある．このような ケースでは，パラメータの値によってどちらのガードが先に満たされるかが異なるため，
離散変化時刻の比較を記号計算によって行う必要が生じるがこれは困難である場合が多 い．ただし非線形モードに限っては，常に非線形関数を包含した平行体上に留まることが できれば，平行体同士の間で遷移を起こすタイミングは自由に選べるはずである．つま り，平行体同士が互いにある程度オーバーラップするような形で包含を行い，平行体間の
第6章 まとめと今後の課題 22 遷移は先にパラメータを分割してから平行体同士がオーバーラップした区間の中で都合の 良い離散変化時刻を選ぶことで，離散変化時刻の比較に伴う問題や離散変化時刻計算によ る誤差の拡大を回避することができると考えられる．
6.2.2 数値シミュレーションを用いた包含の最適化
現在の実装では，単振り子の例題のように場所によって求められる精度が異なるような 場合の考慮や，2変数の微分方程式について軌道が実際には通らないような範囲は重視し ないといった判断を行うことができない．その理由はこれらは実際に微分方程式を解いて みないと分からない情報に依存しており，現在は微分方程式を構成する関数のみを見て包 含の最適化を行っているためである．しかし，実際の運動の予測をもとにして最適化を 行った方が精度は確実に良くなるはずである．そしてこの予測はそれほど厳密に行う必要 はないので，数値計算を用いて事前にシミュレーションを行い，その軌道の速度や通る場 所の情報を利用して最適化を行えばより良い結果が得られるはずである．
6.2.3 三変数以上の微分方程式への対応
現在は線形に変換した微分方程式について，Mathematica を用いて記号解の計算を 行っている．本研究では入力にパラメータを持つことが前提であるため，必然的に微分 方程式も数値ではなく記号のままで解かなければならない部分が多い．これは現在の
Matiematica で行うのはかなり難しく，式によっては 2変数の微分方程式でも解が求め
られない場合が多い．したがって，本手法を単純に三変数以上の場合に向けて拡張しただ けでは，解を求めるのは難しいと考える．
23
謝辞
本研究を進めるにあたり多くのご指導，助言をいただいた上田和紀教授に深く感謝致し ます．また，研究の手法に関して何度も相談に乗っていただき，助言をくださった松本翔 太さんに深く感謝致します．そして，研究室生活の中で議論や相談に乗ってくれた研究室 の皆様に深く感謝致します．
2019年1月 増田 健太
24
参考文献
[1] 上田和紀, 石井大輔, 細部博史: ハイブリッド制約言語HydLa の宣言的意味論，コン ピュータソフトウェア, Vol. 28 No. 1, 2011, pp. 306-311.
[2] Shota Matsumoto and Kazunori Ueda: Symbolic Simulation of Parametrized Hybrid Systems with Affine Arithmetic. In Proc. 23rd International Symposium on Temporal Representation and Reasoning (TIME 2016), Copenhagen, IEEE Computer Society, Oct. 2016, pp. 4-11.
[3] M. Neher, K. R. Jackson, and N. S. Nedialkov: On Taylor Model Based Integra- tion of ODEs, In SIAM J. Numer. Anal., 45 (2007), pp. 236-262.
[4] K. Makino and M. Berz: Suppression of the wrapping effect by Taylor model- based validated integrators, In MSU Report MSUHEP 40910, Michigan State University, 2004.
[5] X. Chen, E. Abraham, and S. Sankaranarayanan: Taylor Model Flowpipe Con- struction for Non-linear Hybrid Systems, In Proceedings of the 2012 IEEE 33rd Real-Time Systems Symposium, pp. 183-192.
[6] Xin Chen, Erika Abraham, and Sriram Sankaranarayanan: FLOW*: An An- alyzer for Non-linear Hybrid Systems, In Computer Aided Verification (CAV), Volume 8044 of LNCS (2013), pp. 258-263.
[7] 柏 木 雅 英: kv ラ イ ブ ラ リ, 早 稲 田 大 学 基 幹 理 工 学 部 応 用 数 理 学 科, 2019, http://verifiedby.me/kv/.
[8] 柏木雅英: ベキ級数演算について, 早稲田大学基幹理工学部応用数理学科, 2019, http://verifiedby.me/kv/psa/psa.pdf.
25
発表論文
[1] 増田健太, 松本翔太, 上田和紀: ハイブリッドシステム解析器 HyLaGI における非線 形常微分方程式の精度保証計算, 情報処理学会第80回全国大会, 2018.
[2] 増田健太, 上田和紀: ハイブリッド制約言語HydLa における非線形常微分方程式の表 現とその記号付き精度保証計算, 2018年度 人工知能学会全国大会（第32回）, 2018.
